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ABSTRACT 


Three  computer  programs  were  written  to  find  the  eigenvalues 
and  eigenfunctions  of  acoustic  normal  modes  in  the  ocean.       The 
programs  used  two  different  methods:      an  iterative  finite  difference 
scheme,    and  a  method  based  upon  the  WKB  approximation  of  quantum 
mechanics.     The  raethods  assume  a  flat  fluid  bottom  and  are  designed 
for  any  arbitrary  sound  speed  profile.      While  the  results  of  both  the 
finite  difference  and  the  WKB  methods  agreed,    the  WKB  method 
proved  faster. 
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I.     INTRODUCTION 

Since  the  Second  World  War  considerable  interest  has  developed, 
in  both  naval  and  scientific  communities,    in  the  prediction  of  sound 
propagation  within  the  ocean.     For  the  scientist,    propagation  informa- 
tion can  be  a  vital  part  in  investigating  the  acoustic,    physical  and 
sometinaes  chemical  properties   of  the  ocean.     For  the  naval  commu- 
nity,   propagation  is  a  vital  element  in  the  prediction  of  acoustic  sen- 
sor performance.     The  prediction  of  such  performance  is  not  limited 
to  the  design  and  development  of  sensor  systems,   but  is  increasingly 
important  in  the  operational  employment  of  such  sensors  and  selec- 
tion of  the  most  fruitful  tactics. 

Accordingly,    there  has  been  considerable  development  within  the 
naval  community  of  numerical  acoustic  propagation  models.      These 
models  have,    for  the  most  part,    been  based  upon  the  theory  of  ray 
acoustics.     The  techniques  employed  have  developed  considerable 
sophistication  in  order  to  deal  with  ray  limitations  with  caustics, 
frequency  effects,    and  interference.     However,    for  prediction  of  prop- 
agation over  great  ranges  and  at  low  frequencies,    ray  techniques 
have  two  serious  limitations: 

(1)    The  long  ranges  involved  require  long  computation  times. 
Each  ray  must  be  "traced"  over  many  successive  short  time  steps  to 
simulate  travel  of  the  wavefront.     Usually  more  than  one  hundred  such 
rays  are  traced.      Thus,    as  ranges  increase,    the  computation  time 
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increases  proportionately.     Such  long  coraputation  times  make  the 
routine  operational  use  of  ray  trace  programs  covering  paths  of 
thousands  of  kilometers  too  expensive  to  be  practical. 

(2)    Ray  trace  programs  largely  ignore  such  frequency  dependent 
effects  as  diffraction  and  interference.      The  assumptions  of  ray  a- 
coustics  require  that  the  wavelength  of  the  acoustic  signal  be  short 
in  comparison  to  the  depth  and  velocity  gradient  scale.     As  lower 
frequency  signals  are  considered,    the  wavelength  becomes  an  appre- 
ciable fraction  of  the  depth  scale,    and  so  violates  a  basic  validity 
assumption  for  ray  acoustics. 

Because  of  these  limitations  there  has  been  increased  interest 
in  normal  mode  theory  as  a  basis  for  long  range  propagation  models. 
The  theory  of  normal  mode  propagation  is  not  new  and  has  had  a  num- 
ber of  shallow  water  applications  [Officer    1958,    Bucker  and  Morris 
1965].     A.    O.    Williams   (1970)  has  outlined  a  method  whereby  the 
normal  mode  technique  could  be  applied  to  deep  ocean  acoustics.     As 
a  part  of  any  such  technique,  both  the  characteristic  horizontal  wave- 
nTimber  (eigenvalue)  and  the  mode  shape  (eigenfunction)  must  be 
solved  for  each  mode.     To  do  this,    three  basic  methods  have  been 
used: 

(1)     The  wave  equation  has  been  solved  in  closed  form  [Tolstoy 
and  Clay   1966,    Williams  and  Home   1967].     In  order  to  do  this,    the 
sound  velocity  profile  must  be  fitted  to  an  assumed  analytic  function. 
Such  functions  have  taken  the  form  of  Epstein  profiles  [Bucker  and 
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Morris   1967],    profiles  with  a  constant  gradient  of  the  reciprocal  of 
the  sound  velocity  or  sound  velocity  squared  [Williams  and  Home 
1967,    Tolstoy  and  Clay   1966].     The  requirement  that  the  sound  veloc- 
ity profile  be  fitted  to  an  arbitrary  function   places  severe  limitations 
on  the  flexibility  of  any  model. 

(2)  The  wave  equation  has  been  vertically  integrated  while  meet- 
ing appropriate  boundary  conditions.      This  method  has  the  advantage 
of  being  adaptable  to  any  arbitrary  sound  velocity  profile.      Kanabis 
(1972)  and  Newman  and  Ingenito   (1972)  have  developed  two  such  shal- 
low w^ater  normal  mode  models.      The  first  program  presented  with 
this  thesis,    NORMOl,    is  based  upon  these  two  shallow  water  models. 

(3)  The  Wentzel-Kramers -Brillouin  (WKB)  approximation  of 
quantum  mechanics  provides  an  analytic  solution  to  the  wave  equation. 
This  solution  is   singular  at  depths  which  correspond  to  ray  vertex 
depths.     However,    the  WKB  approxinaation  may  lend  itself  to  solu- 
tion of  the  eigenvalue,    after  which  the  eigenfunction  may  be  solved 

by  an  integration  similar  to  the  previous  naethod.     The  last  two  pro- 
grams presented  with  this  thesis,    N0RM04  and  N0RM05,    use  this 
method. 

The  long  term  goal  of  this  research  is  to  develop  a  deep  ocean 
normal  mode  acoustic  propagation  model  of  sufficient  computational 
speed  to  be  considered  for  operational  use.     In  order  to  achieve  this 
goal,    a  rapid  method  of  solving  for  the  eigenfunction  must  first  be 
developed.      Based  upon  this  need  the  immediate  objectives  of  this 

thesis  are: 
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(1)  To  compare  the  results  of  the  programs  employing  methods 
(2)  and  (3)  above  with  analytic  solutions,    the  published  results  of 
other  programs,    and  each  other. 

(2)  Determine  whether  the  WKB  method  offers  promise  of  iin- 
provement  in  terms  of  time  versus  accuracy. 

(3)  To  outline  future  improvements  to  either  method  based  upon 
the  preliminary  results. 
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II.     THEORY 

A.      WAVE  THEORY 

In  this  section  we  apply  a  mathennatical  separation  of  variables 
to  the  wave  equation;  then  we  impose  boundary  conditions  on  the  re- 
sulting separated  vertical  equation.      The  results  of  this  raathematical 
procedure  offer  us  a  differential  equation  and  boundary  conditions 
capable  of  numerical  solution.     We  will  close  this   section  with  a  dis- 
cussion of  some  general  properties  of  such  a  solution. 

The  following  discussion  is  based  upon  the  treatments  given  by 
A.    O.    Williams   (1970),    and  I.    Tolstoy  and  C.   S.    Clay  (1966).     For  a 
more  complete  descriptive  treatment  the  reader  is  referred  to  the 
discussion  given  by  Williams. 

1.      The  Wave  Equation 

The  simple  scalar  wave  equation  for  small  amplitude  waves 


IS 


Y72X    _     1    'S'^ 

(1) 


where  <^    is  the  displacement  potential.     Here       ^,    the  displacement 
of  a  fluid  particle  from  its  rest  position,    is  represented  by 


^-  VO, 


(2) 
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and  the  acoustic  pressure,    p,    is  given  by 

Thus,    if  we  can  express  ^    as  a  function  of  space,    we  can  then  find 
the  resulting  acoustic  pressure  and  sound  propagation  loss  due  to 
spreading  and  refraction. 

If  we  use  a  cylindrical  coordinate  system  and  assume  ^  is 
a  function  of  range,    depth,    and  time,    we  can  rewrite   (1)  as 

Let  us  further  assume  that  $   is  separable  in  terms  of  r,    z,    and  t, 
such  that 


^(r,z,t)  -.  R(T-)Z(H)e'"^ 


(5) 


where    <*>    is  the  source  angular  frequency,    Z(z)  is  the  vertical  dis 
placement  potential  function,    and  R(r)  is  the  radial  displacement 
potential  function.     We  immediately  notice  that 


b'§    _    _ 


W^ 


$  .  (6) 


Using  equations   (5)  and  (6)  we  rewrite  equation  (4)  as 


(7) 
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This  expression  has  a  range  dependent  term,    a  depth  depen- 
dent term  and  a  term  with  c,    the  sound  velocity.     If  we  could  assume 
that  c  is  independent  of  either  depth  or  range,    the  expression  is  sep- 
arable into  depth  and  range  differential  equations.      The  obvious   choice 
is  to  assume  that  c  varies  only  with  depth.      This  assumption,    termed 
horizontal  stratification,    although  not  true  over  great  distance,    is 
commonly  made  in  oceanography  and  ocean  acoustics.     Actually,    for 
our  purposes  it  will  be  sufficient  if  c  is  horizontally  stratified  in  the 
vicinity  of  our  solution,    and  further  if  any  horizontal  variation  of  c  is 
small  with  respect  to  the  vertical  variation.     Making  such  an  assump- 
tion,   we  continue  by  separating  (7)  into  range  and  depth  differential 
equations 

2  ^^       F«    "   ^^    •  (9) 

Solutions  to  equation  (8)  include  the  Hankel  functions  of  first 
and  second  kind  representing  cylindrical  waves, 

For  distances  such  that   kj.r    is  much  greater  than  one,    this  function 
can  be  approximated  by  its  asymptotic  form 


--^i^^'^'''--''- 


17 


We  now  have  an  assumed  time  dependence  and  a  solution  to  the  range 
dependent  function  of    SB  .     Thus,    at  long  range  from  the  source  we 
have  for  $ 

J^'kr'^  '  (12) 

2 
Here    kj.    ,    an  eigenvalue  which  appeared  in  the  separation  process  as 

a  mathematical  constant,    will  be  seen  to  be  a  horizontal  wave  number, 

We  w^ill  consider  only  positive  values  of  k    ,    -which  represent  outgoing 

wavefronts.     It  remains  for  us  to  find  a  method  of  solving  for  Z(z)  in 

2 

terms  of  the  separation  parameter,    k      . 

We  can  rewrite  (9)  as 


f5-ff-A?)H>0. 


(13) 


This  equation  is  the  separated,    space  form  of  the  wave  equation,    or 
Helmholtz  equation,    for  the  vertical  wave  component.     The  expression 
in  parentheses  represents  a  vertical  wave  number,    k   ,    and  is  related 
to  the   (true)  wave  number,    k,    and  horizontal  wave  number,    k    ,    by 


r  ■ 


Jk^=^^.Jkl^M^    .  (14) 

Equation  (13)  with  a  pair  of  associated  homogeneous  boundary  condi- 
tions forms  a  Sturm-Liouville  problem,    which  is   solved  in  terms  of  a 
set  of  eigenvalues  and  corresponding  eigenfunctions.     It  now  remains 
for  us  to  define  and  employ  boundary  conditions  associated  with  the 
ocean  vertical  sound  profile.     Before  so  doing,   however,    we  will  di- 
gress a  moment  to  introduce  a  notational  convenience. 
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Since  the  value  of  c  does  not  vary  greatly  (perhaps  five  per- 
cent) throughout  the  depth  profile,    the  square  of  the  vertical  wave 
number  represents  a  small  difference  between  two  nearly  equal  num- 
bers 

Jb^^c^^Jk^,  (15) 

In  order  to  conserve  accuracy  and  facilitate  computational  speed,    a 

notational  convention  is  borrowed  from  quantum  mechanics.      The 

2 
square  of  the  (true)  wave  number,    k    ,    is   subtracted  from  an  arbitrary 

constant  of  the  same  order  of  magnitude   (here  the  maximum  possible 

wave  number)  to  form  a  "potential  function.  " 

YCH)=i#'   -   4'=J.^   -  4'    .  (16) 

The  square  of  the  horizontal  wave  number  (our  separation  constant) 
is  also  subtracted  from  the  maximum  wave  number  to  form  a  quantity 
corresponding  to  the  "energy     level"  of  quantum  mechanics,    and  our 
new  eigenvalue 

E-AL.--tr.  (17) 

With  this  notation,    the  Helmholtz  equation  (13)  becomes 

II    -[e-Y(.)]Z=0   ,  (18) 

a  form  similar  to  Schroedinger 's  equation.     By  this  notation  we  hope 

to  facilitate  the  adaptation  of  the  results  of  quantum  mechanics  to  the 

problem  at  hand. 
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2.      Boundary  Conditions 

We  continue  by  developing  appropriate  vertical  boundary 
conditions. 

Figure   1  represents  a  typical  deep  ocean  sound  profile. 
Figure  2  is  the  same  profile  represented  in  terms  of  V(z).     Note  that 
at  some  intermediate  depth  the  value  of  sound  velocity  reaches  a  min- 
imum at  the  sound  channel  axis.     Note  also  that  the  function  V(z) 
forms  a  "potential  well"  with  depth,    representing  the  deep  sound 
channel.     Our  domain  of  interest  is  bounded  by  the  sea  surface  and  a 
flat  bottom. 

The  air -water  boundary  at  the  surface  approximates  a  free 
surface  at  which  the  stresses  vanish.     In  a  fluid  this  boundary  corre- 
sponds to 

Z(o)  =  0.  (19) 

If  the  limit  of  the  second  derivative  of  the  vertical  displacement  func- 
tion is  finite  at  the  surface,    the  stipulation  that  V(z)  dis continuously 
approaches  infinity  at  the  surface  has  the  effect  of  making  Z{z)  vanish 
at  the  surface.     Thus,    in  terms  of  our  quantum  mechanics  analogy, 
the  free  surface  boundary  condition  corresponds  to  a  perfectly  rigid 
wall  of  the  potential  well. 

The  ocean  bottom  provides  a  boundary  more  difficult  to 
specify.     The  ocean  bottom  has  a  complex,    often  layered,    horizontal- 
ly varying  structure.     In  addition,    little  is  known  of  the  specific 
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TYPICAL  DEEP  OCEAN  SOUND  PROFILE 
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TYPICAL  DEEP  OCEAN  VELOCITY  FUNCTION 
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acoustic  properties  of  much  of  the  bottom  required  tomoclel  fully  the 
bottom  boundary  conditions.     In  the  face  of  such  complexities  we  will 
approximate  the  bottom  with  a  relatively  simple  model. 

We  will  assume  that  the  bottom  is  represented  by  a  horao- 
geneous  fluid  of  infinite  depth,    and  constant  velocity.     Across  a  fluid- 
fluid  interface  we  require  that  both  pressure    and  vertical  particle 
motion  be  continuous.     From  (3)  and  (6),    continuity  of  pressure  re- 
quires 

e.Zo  =  ^bZb  (20) 

at  the  bottom.     From  (2)  we  find  that  continuity  of  vertical  particle 
motion  or  displacement  at  the  bottom  requires 


We  combine  equations   (20)  and(21)  to  form  a  single  bottom  boundary 
condition 

foZo  ^-t        ebZb   a*      *  (22) 

Since  we  have  specified  the  bottom  to  be  of  constant  velocity 
and  density,    a  known  solution  to  equation  (18)  for  the  bottom  region 
is 


Zb-  c.e^«^   ^  c.e-^«^  ,  (23) 


where  Ko,    a  real  constant,    is  the  imaginary  bottom  wave  number, 
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Kb  =i  Vb    -   E  ■    .  (24) 

If  the  amount  of  energy  represented  by  the  bottom  solution  is  to  be 
finite,    the  limit  of  the  integral  of  Z,    from  the  bottom  boundary  to 
infinity  must  in  turn  be  finite.     This  condition  requires  that  the  expo 
nential  growth  term  of  Zi    be  suppressed.      Thus,    (23)  becomes 

b  =  Cze  ,  (25) 

and  (22)  becomes 


ZoTI-'W^B"    ^     ♦  (26) 


We  have  now  specified  two  boundary  conditions  which  will 
enable  us  to  find  particular  solutions  to  the  vertical  differential  wave 
equation  (18).     Before  proceeding  with  the  technique  of  finding  such 
solutions,    let  us  discuss  some  of  their  properties. 
3.       The  Nature  of  the  Solution 

Before  we  find  the  solutions,    Z(z),    to  the  differential  equa- 
tion (18)  and  boundary  conditions   (equations    19  and  26),    it  will  be 
beneficial  to  consider  the  form  and  properties  of  the  solutions  w^e 
expect. 

a.     In  the  previous  section  we  saw^  that  the  solution  in  the 
bottom  region  is  of  the  form  of  a  decaying  exponential  (equation  25), 
This  solution  requires  that  the  quantity  Kg  in  equation  2  5  must  be 
real.     This  requirement  in  turn  implies  that 
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E<Vb  (27) 


or 


kj.)^.  (28) 

c 
b 


In  order  to  have  propagating  waves,    it  must  also  hold  that 

E>0,  (29) 

or 

kj.<^         .  (30) 

min 

Thus,    we  have  an  upper  and  lov/er  bound  for  E  and  k   , 

V^>E>0  (31) 

and 

^      >    K>^    . 

c      .  c  (32) 

min  b 

b.     Within  these  limits  the  boundary  value  problem  is  so  con- 
strained that  it  can  only  be  solved  in  terms  of  a  finite  number  of 
eigenvalues,    E     ,    and  corresponding  functions   Z      (z).     This  eigen- 
function,    properly  normalized,    forms  what  is  termed  a  normal  mode. 
At  sufficient  range  a  vertical  sound  pressure  profile  can  be  approxi- 
mated by  a  linear  combination  of  such  eigenfvinctions, 

tn'i 
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TYPICAL  MODE  SOLUTION 
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Figure  3 
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c.  At  depths  such  that  E^^  V(z)     the  eigenfunction  Z^^(z)  is 
oscillatory  in  character.     We  can  represent  such  a  function  by 

where  \{z)  and     S  (z)    are  undefined  functions  of  depth.     We  shall 
refer  to  such  depths  as  the  oscillatory  or  "real"  region  (a  reference 
to  the  fact  that  the  vertical  wave  number,    k   ,    is  a  real  number).     At 
depths  such  that  E     <  V(z)  the  eigenfunction  Z      (z)  becomes  quasi- 
exponential  in  character,    with  a  form  we  can  represent  by 

where,    again,    ^  (z)  and  S(z)  are  undefined  functions.     We  shall  refer 
to  such  depths  as  the  exponential  or  "imaginary"  region. 

The  depth  at  which  the  sign  of  [e      -V(z)  changes  is  termed 
a  turning  point.     At  the  turning  point  the  character  of  the  eigen- 
function changes  from  oscillatory  to  quasi -exponential     (Fig.    3). 

d.  Each  mode,  corresponds  to  one  or  more  rays  which  traces 
paths  w^ithin  the  oscillatory  region  between  turning  points,    surface, 
or  bottom  boundaries.     The  angle  the  ray  makes  with  the  horizontal 
is  defined  by 

Cos    (Z^Cz)    =   CCz)  4^     .  (36) 

/  CO 

From  this  we  see  that  the  mode  turning  point  corresponds  to  the  ray 
vertex  depth  (that  depth  at  which  a  ray  reaches  its  nnaximum  depth 
excursion).     When  the  oscillatory  region  of  a  mode  is  bounded  by  the 
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free  surface  boundary  or  the  bottom,    the  corresponding  ray  is  sur- 
face or  bottoin  reflected,    respectively. 

e.      The  lower  boundary  placed  on  kj.  represented  by  equation  (28) 
has  a  more  farailiar  and  perhaps  more  satisfying  physical  meaning. 
Substituting  equation  (36)  into  (28)  we  have 


Cos     d   >     C(£)    , 


(37) 


At  the  bottom  this  is  the  expression  for  the  critical  angle.     Thus,    we 
are  limiting  outselves  to  consider  only  those  modes  whose  corres- 
ponding ray  strikes  the  bottom  at  a  grazing  angle  less  than  the  criti- 
cal angle.     Such  modes  are  widely  called  "unattenuated  modes." 
Also,  there  exists  an  infinite  set  of  modes  such  that 


"^      T^  (38) 


However,   because  of  bottom  reflection  loss,    those  modes  tend  to 
attenuate  rapidly  with  range.     For  propagation  problems  at  a  con- 
siderable range,    the  "attenuated  modes"  contribute  an  insignificant 
amount  to  the  acoustic  pressure  field,    and  are  ignored. 

f.     Given  the  oscillatoiry  nature  of  the  eigenfunction  Z     ,    we  see 
that  the  eigenfunctions  for  each  successive  mode  must  change  sign 
one  additional  time  (Figure  4).     Each  sign  change  will  be  referred 
to  as  a  mode  crossing.     Thus,    corresponding  to  each  eigenvalue  E     , 
there  corresponds  an  eigenfunction  Z       with  m-1  mode  crossings. 
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SUCCESSIVE  MODES 


MODE    1 


MODE  2 


MODE  3 


0  crossings 


2  crossings 
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4.     Mode  Parameters 

We  should  now  introduce  two  sets  of  parameters,    derived 
elsewhere,    which  are  used  in  norraal  mode  calculations. 
a.     Phase  and  Group  Speed 

Phase  speed  describes  the  horizontal  speed  of  the  wave 
front  represented  by  a  mode.     As  the  speed  of  advance  of  a  "wave- 
front  of  constant  phase,    the  phase  speed  is  given  by 

Cp  =  ^     •  (39) 

Group  speed  is  the  rate  of  energy  transport  in  the  hori- 
zontal direction.      Tolstoy  and  Clay  (1966)  give  an  expression  for 
group  speed  based  upon  a  theorem  by  Biot  (1957).      This  method  cal- 
culates group  speed  by 


Cg=i.^  (40) 


where  V     is  the  normalization  integral. 


(41) 


and 

-d       ccif  (42) 

b.       Bottom  Decay  Coefficient 

Kornhauser  and  Raney  (1955)  give  an  expression  which 

calculates  the  effect  of  a  bottom  absorption  on  the  mode  amplitude. 

Assume  the  wave  number  in  the  bottom  is   complex  and  given  by 

30 


K-%    *^e  (43) 

where  /9   is  a  bottom  attenuation  coefficient.      Then  the  horizontal 
wave  number  must  also  be  complex  and  given  by 


Ar  =    hv    +  i.  c^  . 


(44) 


If  &   is  small  with  respect  to   —  ,    then  the  ratio  of   o(.    to    (^    for  a  given 
mode  is  approximated  by 


% 


f-^|^/z'<.>d..  (45) 


i.-  boLbom 


B.      THE  WENTZEL-KRAMERS-BRILLOUIN  (WKB)  APPROXIMATION 

The  WKB  approximation  is  a  method  borrowed  from  quantum 
mechanics  which  will  enable  us  to  approximate  the  eigenvalues,    E     , 
of  equation  (18).     In  the  WKB  approximation  we  assume  a  form  for 
the  mode  solution  which  is  valid  at  depths  removed  from  the  turning 
points.     We  then  develop  a  characteristic  equation  for  the  eigenvalue, 
based  upon  the  assumed  solution.      This  section  follows  in  many  re- 
spects the  descriptions  given  by  Tolstoy  and  Clay  (1966)  and  Lauvstad 
(1971). 

1.     The  Characteristic  Equation 

Based  upon  the  sign  of  [E  -  V(z)]  the  vertical  wave  equation 
(18)  applies  to  two  domains,    which  we  have  termed  the  "imaginary" 
region  and  "real"  region.     Let  us  rephrase  equation  (18)  into  two 

separate  equations  for  the  two  domains.     So  doing,    we  have 
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z;  -^  )^i  z.  =  0 


and 


for        E   >   V(E),  (46) 


Z^     -K/   2,   =0 


for        E:   <    VcE).  (47) 


Here  we  have  defined 


K-  [e  -  V(.>] 


Vz 


and 


when       E  >  V(z);,  (48) 


Kz  -    [Yu)  -  E]' 


when       E   <  \Jm.  (49) 


Now  let  us  assume  a  form  for  the  solution  to  equation  (46), 
similar  to  that  which  we  have  previously  postulated  in  equation  (34). 
The  assumed  solution  is 


Zu)  =  5cH)  e 


.5(2) 


(50) 


By  substituting  equation  (50)  into  (46)  we  have  as  a  real  part 


S'%)  -5ci.)C5'fiE)^-J,/)=0, 


(51) 


and  an  imaginary  part 


3cz)S'^Ce)   +  25'z>S'2)   =  0  . 


(52) 


This  equation  immediately  yields  an  expression  for  ^  (z), 

5(H)  -  A/SV/  . 


(53) 
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Now  let  us  simplify  equation  (51)  with  an  assumption.     Assume  that 
the  first  term  of  equation  (51)  is  negligible  with  respect  to  the  others, 


or 


K 


5" 


fz) 


6  w) 


«1 


(54) 


With  this  assumption  we  obtain  as  an  expression  for  S(z) 


S'cH)"  =  Jk^  , 


or  upon  integrating 


il(E)      =     tj 


z 


We  now  have  as  the  assumed  solution 


Zi  (z)  =  X^^  [a-  Cos  Si  (H)   +  b  Sin  Si{e)J, 


(55) 


(56) 


(57) 


or,    in  a  more  convenient  form, 


HifH)  =  i,-'''^[aSLn(5.(H)  +  Q,)] 


(58) 


Similarly  we  obtain  as  a  corresponding  solution  for  the 
"imaginary"  region, 


Z^u,-E;''''[C^^'^'^De-''^'^]^ 


(59) 


where, Sp(z)  is  defined  by 

S2 CH)   =y  Kzda:  . 


(60) 
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Having  ■written  the  two  solutions,    let  us  take  a  closer  look 
at  the  assumption  implied  by  equation  (54).      With  equations   (53)  and 
(55)  we  can  rewrite   (54)  for  the  "real"  region  as 


1_ 


«^    »  (61) 


and  for  the  "imaginary"    region, 


_1_ 


ir''''''\«''  (6z) 


This  assumption  requires  that  the  order  of  magnitude  of  the  vertical 
wave  number,    and  thus  the  sound  speed,    vary  slowly  with  respect  to 
depth.      This  is  generally  the  case  in  the  ocean,    where  the  total  sound 
speed  variation  with  depth  is  on  the  order  of  five  percent.      We  also 
see  that  our  solutions  lose  any  validity  as  z  approaches  a  turning 
point.     Recalling  the  analogy  of  the  turning  point  to  the  ray  vertex 
depth,    we  see  that  as  a  wave  approaches  the  turning  point  its  orienta- 
tion becomes  horizontal.     Near  a  turning  point  the  vertical  wavenum- 
bers,    k^  and  K^,    approach  zero  and  the  conditions  of  equations   (61) 
and  (62)  are  no  longer  satisfied.     Thus  the  solutions  represented  by 
equations   (58)  and  (59)  become  singular  at  the  turning  point  location. 

The  WKB  method  does  not,    at  this  point,    appear  satisfactory 
for  the  computation  of  the  mode  profile  or  eigenfvmction  Z      .     How- 
ever,   if  the  WKB  approximation  can  be  used  to  obtain  an  accurate 
estimate  of  the  eigenvalues,    E     ,    a  finite  difference  scheme  can  then 
be  used  to  calculate  the  eigenfunction,    Z^^. 
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Consider  a  typical  V(z)  profile,    with  upper  and  lower  turning 
points,    as  well  as  bottom  and  pressure  release  boundaries   (Figure  5). 
Remember  that  we  have  specified  a  derivative  boundary  condition  at 
the  bottom.     This  boundary  condition  determines  the  eigenfunction 
(except  for  a  constant  of  multiplication)  in  the  "imaginary"  region 
^TL'^'^b*     "^^  ^^®  lower  turning  point,    ztXj,    the  "imaginary"  region 
solution,    Z2(z),    corresponds  to  a  particular  "real"  region   solution, 
Z,  (z).      We  can  designate  such  a  solution  by  the  use  of  an  initial  phase 
angle,  9  l,    such  that 


(63) 


where 


S^u^-JK^dz  ,  (64) 


BOTTOM 


and 


Sjczj  =  )  .K<^'^-  .  (65) 


z 
2rL 


In  a  similar  manner  we  can  stipulate  that  in  order  to  satisfy  the  free 
surface  boundary  condition,    Z(0)  =  0,    the  "real"  region  solution,    Zj, 
must  correspond  to  a  particular  "imaginary"  region  solution,    Z2. 
Again  we  can  represent  this  as  a  phase  angle,    9g,    in  the  "real"  region 
solution,    Zj.     Thus  we  have  required  the  argument  of  Sin  (Si  (z)  +9j_^) 
at  the  upper  turning  point  (or  surface)  to  be  our  specified  value 
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SiCHru)  *Ql  =■  J  M^c^z  -  Qu  =  Qs   .  (66) 

Ztl 

If  we  let  Q-(x  -"^    "  ®s  ^®  ^^^  rewrite  this  equation  as 

/J^jcIe -G,-0u=  \[l-\/u)]%*Q^*Qu.       =ir.  (67) 

Since  this  equation  describes  the  arguments  for  which  the  sine  of  the 
integrated  vertical  wave  number  is  zero  (a  Dirichlet  boundary  condi- 
tion),   the  surface  boundary  condition  will  also  be  satisfied  if 


r[E-V(H)]''^dE  +  Gl-^Gu   =itiTr  W-1,2,...,N.  (68) 

Ztl 

Values  of  E  satisfying  this  characteristic  equation  will  be  the  eigen- 
values,   E      ,    for  which  we  are  searching.     Note  that  this  equation  is 
also  that  given  by  Tolstoy  and  Clay  (1966)  as  the  characteristic  equa- 
tion for  stratified  acoustic  waveguides. 

Now  that  we  have  a  characteristic  equation,    some  method  of 
evaluating  Q^  and  0     is  required.      The  WKB  approximation  provides 

a  method  for  finding  the  eigenvalues,    E     ,    once  0     and  0,    are  evalu- 

o  o  m  u  L 

ated. 

2.     Turning  Point  Connection 

The  most  vexing  problem  in  using  the  WKB  approximation  is 

connecting  the  two  solutions   (58)  and  (59)  across  the  turning  points.     A 

number  of  techniques  have  been  developed,    and  here  we  will  consider 

two  such  techniques  which  we  will  term  the  asymptotic  and  Airy  phase 

methods . 
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a.     Asymptotic  Method 

We  have  two  corresponding  solutions   represented  by 
equations   (58)  and  (59)  which  correspond  to  the  same  function  on  either 
side  of  a  turning  point.     As  the  turning  point  is  approached  let  us  re- 
quire that  the  limit  of  the  first  derivative  divided  by  the  mode  value 
be  continuous.      This  is  in  fact  the  same  as  the  fluid-fluid  boundary 
condition  of  equation  (22).     For  the  purpose  of  these  turning  point 
discussions  we  will  designate  the  turning  point  depth,    z—     ,    equal  to 
zero  and  corresponding  to  z„  of  equations   (56)  and  (60).     Further, 
z  will  be  positive  towards  the  "real"  region  and  negative  towards  the 
"imaginary"  region  (Figure  6).      Thus,    from  equation  (22)  we  have 

JLim    Zt       =    ilm     Zg    .  (69) 

By  taking  advantage  of  the  approximation  involved  in  equation  (61),    we 
get  by  substitution  for  the  above 


E-Ov    ^  2-0-        Ce*^''^-De-^^'*^  (70) 


Taking  the  limit  we  have 

b.=  CohCGo)  =  D-C.  (71) 

It  should  be  noted  that  this  approach  is  simpler  than  the  Airy  phase 
technique   (described  next)  and  may  not  be  as  accurate.     It  is  included, 


38 


THE  TURNING  POINT  CONNECTION 
A 


Figure  6 
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however,    since  it  provides  an  interesting  comparison  with  the  Airy 
phase  method  in  the  computer  programs. 
b.     The  Airy  Phase  Method 

What  follows  is  an  outline  of  a  technique  developed  by 
V.    R.    Lauvstad  (1971)  of  the  SACLANT  ASW  Research  Centre. 

Lauvstad  assumes  asymptotic  solutions  of  the  same 
form  as  equations   (57)  and  (59)  in  regions   remote  from  the  turning 
points.      The  gradient  of  the  square  of  the  vertical  wave  number  a- 
cross  the  turning  point  is  assumed  to  be  linear,    such  that 

[J^J'=  C.2-....  (72) 

This  reduces  the  vertical  Helmholtz  equation  to 

One  solution  to  this  equation  is  the  Airy  phase  function 

H    =  AAc(-S(z))    -  BB.(-5cE)),  ("^^^ 

where  S(z),    a  depth  integration  function,    is  positive  in  the  "real" 
region,    and  negative  in  the  "imaginary"  region.      The  asymptotic 
forms  of  the  Airy  phase  are  then  compared  with  our  corresponding 
assumed  solutions.     We  can  represent  the  asymptotic  form  of  the 
Airy  pl\ase  for  negative  real  argument  (the  "real"  region)  as 

Zi  =  ±  [(A-B)SinSu)  ^  (A  +  B)Cc.sS(£)]  ,  (75) 
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and  for  positive  real  argument  as 

By  matching  the  above  coefficients  with  those  in  our  solutions  we  have 
a  set  of  coefficient  matching  equations  representing  the  required  turn- 
ing point  connection. 


C=^(oL-b)  D=^(a.b).  (78) 
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III.      THE  PROGRAM  METHODS 

Here  we  will  discuss  the  basic  methods  used  by  the  three  pro- 
grams  (NORMOl,    N0RM04,    and  N0RM05)  presented  as  a  part  of 
this  thesis.      This   section  is  intended  as  only  a  brief  introduction  to 
the  programs  and  their  characteristics.     A  following  section  will  give 
an  outline  of  the  program  algorithms. 

A.      NORMOi    -  A  FINITE  DIFFERENCE  TECHNIQUE 

NORMOl  is  essentially  an  adaptation  of  two  shallow  water  pro- 
grams developed  by  Kanabis   (1972)  and  Newman  and  Ingenito  (1972). 
For  a  given  value  of  the  horizontal  wave  number,    k^^or  E,    the  pro- 
gram begins  with  the  bottom  boundary  condition  (equation  22),    and 
steps  through  a  finite  difference  scheme  to  the  surface.     The  finite 
difference  scheme,    which  is  derived  in  appendix  A,    is  from  Kanabis 
(1972).      The  finite  difference  scheme  is  a  third-order  Newton  forward' 
difference  scheme  which  iterates  both  the  eigenfunction  and  its  deriv- 
ative.    Using  this  difference  scherae,    a  search  is  made  for  those 
values  of  k    ,    or  E,    which  most  closely  satisfy  the  surface  boundary 
condition  (equation  19).     Such  values  are  the  sought  eigenvalues. 

A  shortcoming  of  this  technique  is  that  under  certain  circum- 
stances the  eigenfunction  at  the  surface,    Z      (0),    tends  to  grow  expo- 
nentially rather  than  to  decay  towards  zero.     Some  mention  of  this 
phenomenon  is  made  by  both  Kanabis   (1972)  and  Newman  and  Igenito 
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(1972).     Kanabis   recommends  decreasing  the  incremental  search 
limits  for  the  eigenvalue  (which  imposes  a  consequential  increase  in 
computation  time)  to  overcome  the  phenomenon.     Newman  and  Igenito 
set  the  mode  value  equal  to  zero  at  depths  above  the  last  mode  cross- 
ing or  minimum.     It  happens  that  the  circumstances  for  this  expo- 
nential growth,    or  degenerate  solution,    occur  frequently  in  deep 
ocean  acoustic  profiles.     Thus,    a  more  detailed  treatment  of  the 
phenomenon  is  required. 

To  understand  this  degenerate  solution  let  us  consider  a  typical 
deep  ocean  sound  profile  and  our  general  forms  for  the  eigenfunction 
(equations   34  and  35).      Consider  a  mode  whose  eigenvalue,    E     ,    inter- 
sects the  velocity  function,    V(z),    at  some  depth,    Zj'tt,    below  the  sur- 
face.     This  is  a  common  (in  fact,    the  usual)  occurrence  for  the  lower 
modes  in  the  deep  ocean,    where  a  deep  sound  channel  and  thermocline 
usually  exist.     As  noted  by  Schiff  (1955)  and  Lauvstad  (1971),    if  the 

value  of  S(z)  in  equation  (34)  is  not  exactly -ff^    at  the  turning  point,    the 

4 
coefficients   C  and  D  of  equation  (35)  are  both  non-zero.      Thus,    the 

exponential  growth  term  has   some  finite  value.      Given  both  sufficient 
depth  between  the  turning  point  and  the  surface,    plus  sufficiently 
large  values  of  [V(z)-E],    the  growth  exponential  will  eventually  be- 
come dominant.      Thus,    the  solution  tends  to  grow  exponentially  as  the 
surface  is  approached. 

This  degenerate  solution  can  complicate  the  search  for  eigen- 
values if  it  is  based  solely  upon  the  value  of  Z(0).     With  the 
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degenerate  solution,  the  value  of  Z(0)  will  oscillate  between  extreme- 
ly large  positive  and  negative  numbers,  even  with  only  small  changes 
in  the  horizontal  wave  number  or  E.     For  the  lowest  modes  of  a  deep 

sound  channel  profile  the  values  of  Z(0)  can  so  oscillate  with  incre- 

-14 
mental  changes  in  E  of  only  one  part  in  10         .     Thus,    a  method  such 

as  regula  falsi   (described  later)  can  become  inappropriate  since  a 
map  of  Z(0)  as  a  function  of  E  can  appear  discontinuous   (Figure  7). 
In  order  to  circumvent  this  difficulty,    NORMOl  uses  a  halving 
procedure  based  upon  the  number  of  mode  crossings.     When  the  sur- 
face boundary  condition  is  met,    the  eigenvalue  E  is  the  largest 

7  '  to  j-^  to 

possible  value  of  E  such  that  the  eigenfunction  has  m-1  mode  cross- 
ings.     Based  on  this  property,    the  number  of  mode  crossings  is  used 
to  converge  upon  the  mode  eigenvalue.      The  method  of  successive 
halving,    although  not  sophisticated,    has  the  versatility  to  deal  suc- 
cesfully  with  the  degenerate  solution  and  any  arbitrary  profile.     Once 
an  eigenvalue  is  found,    the  function  Z       is  then  checked  for  a  degen- 
erate solution  near  the  surface.     If  it  is  present  a  new  profile  is 
iterated  in  the  upper  "imaginary"  region  and  matched  to  the  lower 
profile.      This  procedure  may  cause  a  discontinuity  in  the  derivative 
of  Zjj^  at  the  turning  point;  however,    its  ill  effects  would  usually  be 
less  serious  than  those  of  the  degenerate  solution. 

B.   N0RM04  AND  N0RM05  -  THE  WKB  APPROACH 

In  both  N0RM04  and  N0RM05  the  vertical  wave  number  is  inte- 
grated between  the  upper  turning  point  (or  surface)  and  lower  turning 
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point   (or  bottom).      The  programs  differ  only  in  the  turning  point 
connection  formulae.     Both  programs  integrate  the  vertical  wave 
number  using  the  trapezoidal  rule  over  the  depth  grid  points.      The 
eigenvalue  E       is  converged  upon  using  the  regula  falsi  method. 
This  method  makes  successive  estimates  of  the  zero  value  of  a  func- 
tion by  the  linear  interpolation  of  a  negative  and  positive  value.      The 
method  is  represented  by 

El*i   =  E-  +   (E^-E-)(-AS-)    ,  (79) 

(AS^  -  AS-j 

where  the  subscripts   +  and   -  correspond  to  the  values  associated 
with  the  last  positive  and  negative  values  of  A  S;  and  AS,    the  differ- 
ence function,    is  defined  by 

/Ztu 
Jk^6z  ^Gu-e.-'m'^-  (80) 

Htv 

If  the  mode  is  not  found  after  a  set  number  of  iterations,   the  regula 
falsi  method  is  dropped  for  a  halving  process. 

The  phase  effect  of  the  upper  and  lower  turning  points,    9^  and 
0,  ,    are  evaluated  using  the  asymptotic  method  for  N0RM04  and  the 
Airy  phase  method  for  N0RM05.     The  vertical  wave  number  is  integ- 
rated from  the  uppermost  turning  point  (or  surface,    if  surface  re- 
flected) to  the  lowest  turning  point  (or  bottom,    if  bottom  reflected). 
Both  programs  must  then  provide,    in  addition  to  the  values  of  0^  and 
Qj,   the    effect  of  any  intermediate  "inaaginary"  region  or  barrier. 
This  occurs  when  multiple  sound  channels  are  encountered  (Figure  8) 
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Details  of  the  derivation  of  the  specific  formulae  for  9    ,    9t  ,    and  the 
phase  effect  of  an  intermediate  barrier  are  given  in  appendices  B  and 
C  for  N0RM04  and  N0RM05  respectively. 

If  the  mode  is  surface  reflected  or  bottom  reflected,   both  pro- 
grams use  the  same  formulae  for  the  phase  effect  of  the  reflection, 
0^  and  0T  .     Since  the  free  surface  requires   Z(0)  =  0,    for  the  surface 
reflected  wave,    we  have  from  equations   (63)  and  (66) 


L 


«<,dz  +6^  =  9^=  tnTf  (81) 


and  thus 


9.-0.  (82) 

For  the  bottom  reflected  mode,    we  require  the  conditions  of  a  fluid 
fluid  boundary  to  be  met  (equation  22).     Substituting  equations   (63) 
and  (25)  into  (22)  we  have 


J,,CotC8.V--^KB   ,  (83) 


or 


where  k^  is  evaluated  at  the  water  side  of  the  bottom  boundary. 

The  occurrence  of  multiple  sound  channels  or  ducts  makes  the 
WKB  method  more  complicated.     If  the  intermediate  "imaginary" 
region  or  barrier  is  of  sufficient  extent  and    magnitude,    sound 
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propagation  in  each  channel  becomes  independent  of  the  other.     When 
this  is  the  case,    the  waves  in  one  channel  undergo  an  effectively  total 
reflection  by  the  barrier.      Thus,    complete  sound  ducting  occurs 
entirely  within  that  channel.     Lauvstad  (1971)  gives  the  ratio  of  the 
amplitudes  of  the  reflected  and  incident  waves  in  a  duct  as 

Ri=   e*-  -  e"*-  (85^ 

e*-  +■  e"^ 

where  L  is  the  integral  of  the  imaginary  vertical  wave  number  across 
the  barrier. 


(86) 


This  result  can  be  derived  using  either  the  asymptotic  or  Airy  phase 
connection  formulae.     As  L  increases,    Ri   rapidly  approaches  unity, 
and  near  total  reflection  is  achieved  by  the  barrier.      The  WKB  pro- 
grams    N0RM04  and  N0RM05     evaluate  L,,    and  if  it  is  of  sufficient 
magnitude,    the  programs  treat  the  mode  propagation  in  each  channel 
independently.     The  programs  accommodate  a  maximum  of  one  main 
sound  channel  (which  contains  the  velocity  minimum),    and  one  upper 
and  one  low^er  sound  channel. 
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IV.      THE  BASIC  ALGORITHMS 

A.      NORMOl 

1 .       General  Description 

NORMOl  iterates  the  wave  equation  vertically  from  t?ie 
bottom  to  the  surface  using  a  finite  difference  scheme  over  a  depth 
grid  of  up  to  501  points.     The  program  is  written  in  IBM  FORTRAN 
IV  with  subroutines  used  as  major  functional  blocks.      This  modular 
programming  is  designed  for  both  ease  of  modification  and  debugging. 
A  description  of  the  sequence  of  events  follows. 

NORMOl  is  basically  composed  of  three  nested  loops:     a 
frequency  loop,    a  mode  loop,    and  an  iteration  loop.      The  input  data 
are  read  using  an  input  subroutine.     The  required  input  data  are 
listed  in  appendix  D.     Once  the  data  are  read  and  appropriate  con- 
stants computed,    the  sound  velocity  profile  is  linearly  interpolated 
for  all  grid  points.     Here  the  frequency  loop  begins,    the  appropriate 
frequency  is  selected,    and  the  velocity  function  V(z)  is  computed  for 
the  grid  points.     Next  the  maximum  and  minimum  horizontal  wave 
numbers  and  maximum  number  of  modes  are  computed.     If  the  fre- 
quency is  below  cut-off  (the  maximum  number  of  modes  is   zero),    a 
new  frequency  is  selected.     At  this  point  the  mode  loop  begins,    and 
the  following  steps  are  repeated  for  the  first  to  the  highest  requested 
or  possible  mode. 
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First,    a  method  of  successive  halving  is  used  to  find  two 
values  of  the  horizontal  wave  number  k     which  have  m-1  and  m  mode 
crossings.      The  desired  eigenvalue  must  lie  betv/een  these  two  values, 
Then,    again  using  a  halving  process,    the  program  converges  upon  the 
eigenvalue  k    .      The     mode  is   considered  found  when  one  of  two  con- 
ditions   (regulated  by  an  input  parameter,    lEX)  is  satisfied; 

IZcoil    i    |Z1„„.10"'",  (87) 

or  the  horizontal  wave  number  changes  between  iterations  by  an 
amount  such  that 


^. 


<    10  -^^^^  .  (88) 


Once  the  mode  is  found,    the  eigenfunction  is   checked  for  a 
degenerate  solution  (exponential  growth  near  the  surface).     If  a  de- 
generate solution  is  present,    a  new  function  is  iterated  from  the  sur- 
face to  the  upper  turning  point,    then  joined  w^ith  the  lower  function. 
The  mode  is  then  normalized  to  its  absolute  maximum  value,    and  the 
profile  within  the  homogenous  bottom  is  computed.     The  computed 
mode  is  plotted  using  a  printer  plot  routine.     A  subroutine  then  com- 
putes the  normalization  integral  represented  by  equation  (41),    and 
the  group  speed,    phase  speed  and  mode  attenuation  coefficient.     At 
this  point  the  mode  and  frequency  loops  end. 

Included  in  NORMOl  is  an  output  subroutine  which  provides 
a  punched  card  format  listing  of  the  eigenvalues,    mode  parameters. 
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sound  velocity  profile  and  naode  profile.     The  output  subroutine  was 
written  to  provide  local  input  to  CALCOMP  plots  and  propagation 
loss  prqgrams. 

2.     The  Subroutines 

The  following  is  a  brief  description  of  the  major  subroutines 
which  comprise  the  NORMOl  program. 

INPUT  1  provides  the  input  data  for  the  program.     A  user 
supplied  page  heading  is  read  on  the  first  card,    followed  by  the  run 
parameters,    list  of  frequencies  desired,    and  a  maximum  of  50  depth 
and  sound  velocity  pairs.     If  the  parameter  representing  the  number 
of  depth  and  velocity  pairs,    NUMV,    is  negative,    all  length  measure- 
ments read  are  converted  from  feet  to  meters.      This   subroutine  also 
prints  the  input  data  for  reference. 

INTRPL  linearly  interpolates  the  input  depth  and  sound 
velocity  pairs  to  compute  the  grid  point  values  of  sound  velocity.     If 
required,    this  subroutine  also  extrapolates  the  last  sound  velocity 
to  the  bottom  assuming  an  isothermal,    isohaline  gradient  (0.017sec~  ). 

MAXMIN  computes  and  checks  the  raaximum  and  minimum 
values  of  the  horizontal  wave  number,  and  finds  the  number  of  modes 
represented  by  the  minimum  horizontal  wave  number  (the  maximum 
number  of  modes).  If  the  maximum  number  of  modes  is  zero,  the 
frequency  is  below  cut-off  and  the  program  selects  a  new  frequency 
or  terminates  the  calculations.  If  the  maximum  possible  number  of 
modes  is  below  that  requested,    this  parameter  is    changed  and  noted. 
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HALF  searches  for  two  values  of  the  horizontal  wave  number 
with  m  (the  mode  being  searched)  and  m-1  mode  crossings.     These 
values  represent  lower  and  upper  bounds  on  k     which  are  then  suc- 
cessively narrowed.     At  the  sarae  time  this   subroutine  "maps"  the 
values  of  k     versus  mode  number  in  order  to  facilitate  the  search 
for  later  modes. 

ITERAT,    the  heart  of  the  program,    iterates  the  eigenfunc- 
tion  from  the  bottom  to  the  surface  using  equations   (97)  and  (101). 
In  addition,    this  subroutine  counts  the  number  of  mode  crossings  and 
finds  the  maximum  absolute  value  of  the  eigenfunction. 

DNORMl  normalizes  the  eigenfunction  with  respect  to  its 
maximum  absolute  value. 

BOTTOM  creates  an  exponentially  decaying  bottom  mode 
profile  from  the  bottom  to  an  input  depth,    PLTMAX. 

MODPLT  and  CZPLOT  are  printer  plot  subroutines  which 
plot  the  mode  shape  and  sound  velocity  profile.      Both  these  subrou- 
tines use  a  Naval  Postgraduate  School  utility  printer  plot  subroutine, 
UTPLOT. 

FIXUP  tests  for  a  degenerate  solution  near  the  surface.     If 
the  degenerate  solution  exists,    this  subroutine  applies  the  finite  dif- 
ference iteration  scheme  starting  from  the  surface  downward  to  the 
upper  turning  point.     This  upper  raode  shape  is  then  scaled  to  join 
the  original  lower  function  at  the  turning  point.     If  required,    a  new 
value  for  the  absolute  maximum  of  the  eigenfunction  is  then  found. 
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INTEGR  computes  the  mode  normalization  integral  (equa- 
tion 41)  and  the  integral  represented  by  (equation  42).      The  phase 
velocity  is   computed;    then  the  two  integrals  are  used  to  compute  the 
group  velocity.     Finally,    the  mode  attenuation  coefficient  (equation 
45)  is  calculated. 

B.      N0RM04  AND  N0RM05 
1.     General  Description 

N0RM04  and  N0RM05  are  identical  with  the  exception  of  the 
turning  point  connection  formulae  (appendices  B  and  C  explain  these 
differences).     Similar  to  NORMOl,    they  consist  of  three  nested  loops; 
a  frequency  loop,    a  mode  loop,    and  an  iterative  loop  represented  by 
the  subroutine  SEARCH.     In  addition,    the  input,    interpolation,    plot- 
ting,   maximum /minimum  and  mode  integration  functions  are  per- 
formed by  subroutines  similar  to  those  used  in  NORMOl. 

The  sequence  of  events  in  the  WKB  method  programs  is 
similar  to  NORMOl.     N0RM04  and  N0RM05  read  the  input  data  and 
then  interpolate  the  sound  velocity  profile  over  a  maximum  of  1001 
depth  grid  points.     If  requested,    a  sound  velocity  printer  plot  or 
punched  card  format  sound  profile  is  provided.     The  frequency  loop 
then  computes  the  velocity  function  V(z)  and  checks  for  maximum 
horizontal  wave  number,    minimum  horizontal  wave  number,    fre- 
quency cut-off,    and  the  maximum  number  of  modes  available.     At 
this  point  the  mode  loop  begins  and  the  search  for  the  eigenvalue  is 
controlled  by  the  SEARCH  subroutine. 
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After  the  mode  is  found,    the  mode  shape  is  computed  using 
a  finite  difference  scheme.     The  finite  difference  scheme  is  fitted  be- 
tween the  upper  and  lower  turning  points,    then  a  solution  which  decays 
away  from  the  turning  points  is  calculated.     After  the  mode  shape  has 
been  calculated  and  normalized,    the  mode  is  integrated  to  provide  the 
normalization  integral,    phase  and  group  speed,    and  the  mode  attenu- 
ation coefficient.     Finally,    if  requested,    a  printer  plot  of  the  mode 
shape  is  provided.     At  this  point  the  mode  and  frequency  loops  end. 

Because  of  their  importance  to  these  programs,  the  SEARCH 
and  WKB  subroutines  deserve  more  extensive  description. 
2.     SEARCH  Subroutine 

The  SEARCH     subroutine  controls  the  iterative  search  for  the 
eigenvalues  k    .      The  subroutine  first  makes  an  initial  guess  of  an 
upper  and  lower  bound  for  the  eigenvalue.      The  WKB  subroutine   is 
then  called  to  compute  the  difference  function  given  by  equation  (80). 
If  the  difference  function  for  the  lower  bound  is  not  positive,    the 
bound  is  successively  relaxed  until  a  positive  value  is  found. 

The  next  guess  for  the  eigenvalue  is  calculated  using  a 
regula  falsi    zero -finding  technique  (given  by  equation  79).     The  re- 
sulting incremental  change  in  eigenvalue  and  the  absolute  value  of  the 
difference  function  is  checked  to  test  whether  the  desired  accuracy 
has  been  achieved.     The  WKB  subroutine  is  then  called  to  find  the 
difference  function  for  the  new  trial  eigenvalue.     If  the  resulting  dif- 
ference function  is  negative  the  upper  bound  is  replaced;  if  it  is 
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positive  the  lower  bound  is  replaced.     A  new  guess  for  k     is  then 
estimated,    and  this  iterative  process  continues  until  one  of  two  con- 
ditions is  met.      The  mode  is  considered  found  when  the  difference 
function  falls  below  an  absolute  limit, 

|AS|     <    10  "^^^  '         (89) 

or  the  value  of  the  horizontal  wave  number  approaches  machine 
accuracy, 


<   10 


-14  (90) 


If  after  a  set  number  of  iterations,    the  regula  falsi  method  has  failed 
to  converge  upon  a  solution,    it  is  replaced  by  a  cruder  and  slower 
(but  more  versatile)  halving  process. 

If  more  than  one  sound  channel  exists,    the  program  must 
determine  if  the  two  channels  are  independent  of  one  another.     If  the 
L  of  equation  (86)  is  of  sufficient  magnitude,    the  subroutine  SEARCH 
computes  the  eigenvalues  for  each  channel  separately.     Thus,    after 
an  eigenvalue  has  been  found,    two  tests  are  made.     First,    it  is  deter- 
mined -whether  an  upper  or  lower  (separated)  sound  channel  exists. 
Then,    if  an  upper  or  lower  channel  exists,    it  is  tested  to  determine 
whether  it  will  propagate  an  additional  mode.     If  such  a  mode  can  be 
propagated,    the  SEARCH  subroutine  computes  its  next  mode  within 
that  channel. 
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3.       WKB  Subroutine 

The  purpose  of  the  WKB  subroutine  is  to  integrate  the  verti- 
cal wave  number  between  the  upper  and  lower  turning  points.      The 
subroutine  consists  of  a  loop  which  steps  through  the  depth  grid  from 
the  surface  to  the  bottom.     In  any  "imaginary"  region  the  imaginary 
vertical  wave  number  K     is  integrated  for  use  with  the  turning  point 
connection  formulae.     Each  time  the  depth  loop  exits  an  "imaginary" 
region  for  a  "real"  region,    either  the  upper  turning  point  phase  0 
or  the  phase  effect  of  a  barrier  must  be  calculated.     If  the  barrier  is 
considered  to  be  separated  (if  L  in  equation  86  is  larger  than  a  set 
value)  the  integration  is  stopped  with  the  upper  channel. 

When  the  depth  loop  reaches  the  bottom,    the  bottom  phase, 
Q-r  ,    is  computed  for  either  a  bottom  reflected  or  a  decaying  solution 
as  appropriate.     Finally,    the  difference  function  is  calculated  and 
program  control  returns  to  the  calling  program. 
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V.     RESULTS 

All  three  programs  were  run  with  a  variety  of  sound  velocity 
profiles  and  the  calculated  eigenvalues  examined.      The  profiles  select- 
ed   include  the  published  results  of  other  models  and  profiles  with 
known  analytic  solutions.     In  all  the  test  runs  the  input  variable  lEX 
was  set  at  10.      Thus,    the  iterations  were  halted  when  the  absolute 
relative  value  of  the  eigenfunction  at  the  surface,      Z(o)    ,    or  the  dif- 
ference function,   A  S    ,    (equation  80)  achieved  a  value  less  than  10 
Additionally,    computation  was  halted  if  the  eigenvalue  was  incremented 
less  than   10  during  a  halving  process   (in  other  words,    when  the 

solution  approached  machine  accuracy). 

A.      TEST  CASE  ONE  -  NEGATIVE  GRADIENT 

This  test  case  is  one  given  by  Kanabis   (1972).     The  profile  is 
characterized  by  a  simple  negative  velocity  gradient  in  a  shallow 
water  channel  (Figure  9).     The  profile  has  a  fast  fluid  bottom  of  the 
same  density  as  the  water.     The  test  results  are  computed  for  a  fre- 
quency of  1000  Hz.     Table   1  lists  the  eigenvalues  given  by  Kanabis 
and  computed  by  NORMOl,    N0RM04  and  N0RM05  in  yd      .      The  list 
includes  a  set  computed  by  the  Kanabis  program  and  a  set  computed 
by  an  unpublished  program  by  C.    L.    Bartberger  and  L.    E.    Ackler. 

Except  for  the  higher  nnodes  of  N0RM05,    the  results  from  all  three 

5 
programs  agree  to  the  fourth  decimal  place  or  one  part  in  10    . 
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TEST  CASE  ONE  SOUND  PROFILE 
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Table  2  gives  the  NORMOl,    N0RM04  and  N0RM05  results  to 
the  accuracy  computed  and  relative  differences  between  NORMOl  and 
the  other  two  prograras.     NORMOl  required  56.23  sec  of  central 
processing  unit  (CPU)  time  and  412  iterations  to  arrive  at  the  solu- 
tions.      N0RM04  and  N0RM05  required  23.  50  seconds  with  83  itera- 
tions and  22.  74  seconds  with  79  iterations  respectively.     The  CPU 
times  cited  include  the  time  required  for  all  operations  including  in- 
put,   sound  speed  profile  manipulation  and  graph  printing. 

B.      TEST  CASE  TWO  -  POSITIVE  GRADIENT 

The  second  test  case  is  a  shallow  water  example  cited  by  Newman 
and  Ingenito  (1972).     The  profile  is  characterized  by  a  slow,    variable 
positive  sound  speed  gradient  over  a  fast,    dense  bottom  (Figure   11). 
The  computation  frequency  was  700  Hz.     It  should  be  noted  that  the 
accuracy  parameter  used  by  Newman  and  Ingenito  was  an  absolute 
value  of  less  than  0.001  for  the  surface  eigenfunction  Z(0).     This 
roughly  corresponds  to  a  value  of  lEX  between  3  and  4  in  the  NORMO 
programs.      With  the  exception  of  N0RM05,    all  the  programs  agreed 
to  at  least  the  third  decimal  place  (Table  3).      Considering  the  accu- 
racy limit  set  for  the  Newman  and  Ingenito  run,    this  agreement  is 
about  as  good  as  may  be  expected.     NORMOl  required  27.43  seconds 
and  196  iterations  to  complete  computations,    while  N0RM04  and 
N0RM05  required   14.  73  seconds  for  77  iterations  and   14.  92  seconds 
for  79  iterations,    respectively. 
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TEST  CASE  TWO  PROFILE 
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C.      TEST  CASE  THREE   -  SYMMETRIC  WAVE  DUCT 

A  profile  was  constructed  corresponding  to  a  symmetric  wave 
duct  where  the  sound  speed  is  given  by 


=  JL   _  a^w^fzo-z] 


Z  /-i 


c      c, 


o 


(115) 


between  two  depths.     The  solution  to  this  profile  is  given  by  Tolstoy 
and  Clay  (1966)  as 

""        L  Co*-  J 

The  actual  profile  dimensions,    shown  in  figure   12,    are  taken  from 
Tolstoy  and  Clay.      The  profile  was  tested  at  30.0  and  60.0  hertz. 
The  computed  eigenvalues  agree  with  the  anlytic  solutions  to  one 
part  in  10    ,    with  the  exception  of  the  higher  modes  of  NORMOl. 
Table  4  shows  comparisons  betw^een  the  analytic  solution  and 
NORMOl,    N0RM04  or  N0RM05  computed  values. 

Of  significant  interest  is  the  fact  that  the  errors  appeared  rela- 
tively stable  for  N0RM04  and  N0RM05.     The  errors  for  N0RM04 
are  consistently  0.  45x10       to  0.69x10"     meters  ~      low  for  30  Hz,    and 
0.  lOxlO"^  to  0.  12x10"^  meters  "     low  for  60  Hz.     In  contrast,    the 
errors  for  NORMOl  increased  from  0.  39x10"     meters  "     low  to 
12.  5x10"^  meters  "^   high  for  30  Hz,    and  from  0.  SOxlo"     to  2.  8x10'^ 
meters  ~ '■  low  for  60  Hz.     Figure   14  displays  a  graph  of  the  program 
errors.      Figure    13  shows  the  first  three  modes. 


66 


TEST  CASE  THREE  PROFILE 


f>liG 


— I 


Mcj 


ZT 


1000. 0  m 
1539       m/sec 


m/sec 


Z2=  2828.8m 
Cy-  1539      m/sec 


BOTTOM 


Figure   12 
67 


TEST  CASE  THREE  MODES 
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D.   TEST  CASE  FOUR  -  DEEP  SOUND  CHANNEL 

Test  case  four  is  a  typical  deep  ocean  sound  channel.      The  pro- 
file is  characterized  by  a  deep  and  sharp  sound  channel  and  by  a 

bottom  of  the  saine  density  as  the  water  (Figure   15),      With  the  excep 

4 
tion  of  the  first  raode,    the  eigenvalues  agree  within  one  part  in  10 

(Table  5).     Note  that  for  the  higher  frequency,    60  Hz,    the  agreement 

is  somewhat  improved.     It  is  interesting  that  the  first  mode  should 

show  a  larger  disagreement  among  the  programs.     This  effect  may 

be  a  result  of  the  combination  of  the  sharp  sound  channel  axis  and  a 

grid  spacing  of  8  meters. 
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TEST  CASE  FOUR  PROFILE 
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RESULTS  OF  TEST  CASE  FOUR 


Frequency:     30  hertz 

NORMOl         N0RM04        xlO^ 


MODE 


4 
5 
6 
7 
8 
9 
10 


0.  12614128  0.  12610844 
567736 
537387 
514437 
492671 
472443 
453745 
436162 
419601 
404106 


Frequency:  60  hertz 

1      0.25255354   0.25250819 


4 
5 
6 
7 
8 
9 
10 


196659 


154658 


118755 


089177 


063654 


039561 


016594 


196827 


154783 


118986 


089755 


039343 


016618 


0.24994761   0.24994890 


N0RM05 


DIFF 
xlO^ 


567933 

-  1.  97 

538187 

8.00 

513923 

-  5.  14 

492224 

-  4.47 

472294 

-  1.49 

453659 

-  0.86 

436012 

-  1.50 

419117 

-  4.84 

400638 

-34. 68 

32.84      0.12610844      -32.84 


567933       -    1.79 


538187  8.00 


513923      -  5.  14 


492224      -  4.47 


472294       -   1.49 


453660      -  0.85 


436018 


1.44 


419155      -  4.46 


402169      -19.37 


■45.35      0.25250819      -45.35 


1.  68 


1.25 


2.  31 


5.78 


196827 


154783 


118986 


1.68 


1.25 


2.31 


063571  -  0.83 


2.  18 


0.24 


089755  5.78 


063571      -  0.83 


039343      -  2.  11 


016618  0.24 


1.29      0.24994890  1.29 


973743 


974255 


5.  12 


974255  5.  12 


TABLE  V 
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E.      TEST  CASE  FIVE  -  DOUBLE  CHANNEL 

Test  case  five  is  a  double  sound  channel  of  nearly  symnaetric 
dimensions   (Figure   16).     This  profile  was  used  to  test  the  multi-duct 
properties  of  the  N0RM04  and  N0RM05  programs.     Note  that  for 
30  Hz  the  differences  between  NORMOl  and  N0RM04  (or  N0RM05) 
increase  by  almost  an  order  of  magnitude  for  modes  5  and  7  (Table  6) 
For  these  modes  the  barrier  connection  formula  (equation  114)  is 
employed  and  the  ducts  are  considered  "connected."     The  mode  pro- 
files  (Figure   17)  for  these  modes  also  do  not  completely  correspond. 
Modes  5  and  7  appear  to  be  "upside-down"  in  N0RM04,    with  the 
dominant  amplitude  waves  in  the  lower  channel.     This  is  due  to  the 
fact  that  N0RM04  and  N0RM05  iterate  the  vertical  wave  number 
from  the  surface  downward.     Because  of  this,    the  phase  at  the  upper 
turning  point  of  the  lower  channel  represents  the  phase  angle  for  the 
third  and  fourth  mode  in  the  lower  channel.     Equation  (114)  remains 
correct,    although  here  it  has  been  incorrectly  applied. 
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TEST  CASE  FIVE  PROFILE 


scijNri  (.Ei-aniT"- 


fXiQ 


i.SQ 


BOTTOM 


Figure  16 
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RESULTS  OF  TEST  CASE  FIVE 


MODE 


NORMOl 


N0RM04 


DIFF 
xlO^ 


N0RM05 


Frequency:     30  hertz 

1  0. 12608464      0. 12604450 


2 
3 

4 
5 
6 


08403 


04450 


0.12553486   0.12554660 


52366 


16376 


11170 


53937 


36343 


09716 


0.12485241   0.12490245 


10 


75394 


59913 


42282 


74793 


59806 


42565 


Frequency:  60  hertz 

1     0.25248206   0.25242998 


48204 


42998 


0.25179326   0.25180252 


79260 


31541 


30922 


80252 


31134 


31100 


7 

8 

9 

10 


0.25087788   0.25088532 


84379 
50385 
41349 


84790 


49918 


40770 


-39.53 


04450 


11. 76  0.  12554660 


15.  71 


199. 67 


-14. 54 


53937 


50. 04  0.  12490245 


-  6.  01 


-  1.07 


2.83 


74793 


59806 


42568 


•52.08  0.25242998 


-52.06 


42998 


9.26  0.25180252 


9.92 


-  4.07 


1.78 


80252 


31134 


31100 


7.44  0.25088532 


4.  11 


-  4.  67 


5.  79 


84790 
49918 
40770 


TABLE  VI 
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DIFF 
xlO^ 


■40. 14    0. 12604450         -40. 14 


■39.53 


11.76 


15.  71 


36343         199.67 


09716         -14.54 


50.04 

•  6.01 

■  1.07 
2.86 

■52/08 

■52.06 

9.26 

9.92 

•  4.07 

1.78 
7.44 
4.  11 

■  4.  67 

•  5.  79 
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'igure    17 
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VI.     CONCLUSIONS 

A.      GENERAL 

As  seen  by  table  7,    the  WKB  method  has  performed  consistently- 
faster  than  the  iterative  technique.      The  WKB  programs    (N0RM04 
and  N0RM05)  found  solutions  with  about  one -half  the  Central  Pro- 
cessing Unit  (CPU)  time  required  by  the  finite  difference  program 
(NORMOl).     In  addition,    one  should  consider  the  method  used  for  in- 
tegration of  the  vertical  wave  number  in  the  WKB  programs.      The 
trapezoid  rule  was  used  which  required  the  evaluation  of  a  square 
root  at  each  grid  point.     Although  simple,    this  is  a  time  consuming 
task,    while  it  is  not  necessarily  very  accurate.     If  a     faster  integra- 
ting procedure  w^ere  substituted  for  the  trapezoid  rule,    the  time  re- 
quired for  each  iteration  could  be  reduced  considerably.      Thus,    in 
addition  to  showing  an  improvement  in  computation  time,    the  WKB 
method  also  offers  the  potential  of  greater  computational  speed 
through  the  use  of  sophisticated  integrators. 

The  WKB  and  finite  difference  methods  appear  to  have  similar 
accuracy  in  comparison  with  other  programs.     The  accuracy  attain- 
able is  adequate  to  calculate  useful  modes  for  a  propagation  loss  pro- 
gram.    However,    the  errors  involved  in  the  differences  between 
successive  eigenvalues  will  not  permit  accurate  coherent  phasing  of 
the  modes  beyond  the  first  or  second  convergence  zone    .      Coherent 
phasing  at  extremely  long  ranges  is,    in  fact,    attainable  only  with  an 
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CPU  TIMES   FOR  VARIOUS  RUNS 


Test  Case  One 


Test  Case  Two 


Test  Case  Four 


Test  Case  Five 


NORMOl 


56. 23  sec 


27.43 


Test  Case  Three  106.  35 


105. 67 


114.07 


N0RM04 
23.50  sec 
14.73 
38.42 
40.23 
44.53 


N0RM05 
22. 74  sec 
14.  92 
39.47 
40.43 
47.09 


TOTAL 


409. 75  sec 


161.41  sec 


164. 65  sec 


TABLE  VII 
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analytic  profile  and  solution  —    an  ideal  mathematical  model.     In 
fact,    when  one  considers  the  accuracy  of  the  available  sound  velocity 
data,    it  is  difficult  to  imagine  any  model  attempting  to  yield  accurate 
phasing  effects  at  any  large  distance.     Since  the  position  of  a  target 
is  unknov/n,    in  practice  it  is  more  important  that  the  amplitude  and 
nature  of  phasing  effects  be  know,    than  for  the  sharp  interference 

peaks  to  be  accurately  predicted.      Thus,    accuracies  of  one  part  in 

5 
10     should  enable  a  model  both  to  take  account  of  the  phasing  effects 

in  the  first  convergence  zone,    and  to  give  a  general  description  of  the 
effects  of  phase  interference  at  longer  ranges. 

Although  the  WKB  programs  do  not  correctly  interpret  the  barrier 
effects,    the  problem  appears  to  lie  in  the  approach  taken  in  the  pro- 
gramming,   rather  than  in  the  mathematical  and  physical  relationships. 
Within  any  duct,    the  effects  of  any  upper  and  lower  ducts  must  be 
treated  as  a  turning  point  phase  effect  similar  to  the  top  and  bottom 
boundary  conditions.     The  modes  5  and  7  of  the  N0RM04  30  Hz  double 
channel  case  are  unwitting  examples  of  this.     In  N0RM04  and 
N0RM05  the  "connected"  ducts  are  considered  as  a  single  wave  sys- 
tem.    In  fact,    with  any  sized  barrier,    no  matter  how  small,   one  must 
consider  the  wave  system  within  each  duct  separately.     In  this  case 
the  turning  point  phase  angle  between  the  duct  and  barrier  is  analo- 
gous to  the  impedence  due  to  the  other  duct  felt  by  the  wave  system  in 
its  duct. 
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B.      FUTURE  PROGRAM 

It  appears  that  a  faster  WKB  method  normal  mode  program  can 
be  written  which  will  properly  interpret  a  barrier.     In  order  to  pro- 
perly consider  the  multiple  channel  cases,    a  mapping  procedure 
should  be  used  before  actual  mode  searches   commence.      The  sound 
profile  would  be  divided  into  "levels"   (Figure   18),    found  by  selecting 
all  the  velocity  raaxima  and  minima.     Each  "level"  represents  a 
range  of  the  possible  eigenvalues,    E      .      Within  each  "level"  there  is 
a  single  combination  of  ducts,    barriers  and  boundaries.      At  the  start 
of  computations  the  program  would  determine  the  limits  of  each  level. 
Then,    at  the  beginning  of  the  frequency  loop,    the  program  would  de- 
termine the  maximum  and  minimum  number  of  modes  for  each  duct 
within  each  level.     Thus,    equipped  with  a  map  of  the  profile  charac- 
teristics,   the  program  would  compute  eigenvalues  for  each  duct  and 
level  in  sequence   (Figure   19).     After  all  the  eigenvalues  are  computed 
and  sorted  (some  may  be  out  of  numerical  order),    the  eigenfunctions 
can  be  computed  quickly  using  the  finite  difference  procedure  of 
NORMOl.     Combined  with  a  faster  integrating  method,    such  a  pro- 
gram should  yield  fast,    accurate  normal  modes  for  an  arbitrary 
sound  speed  profile. 
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SOUND  VELOCITY  PROFILE 
WITH 
LEVELS 


velocity   funcfion 


Figure    18 
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SEQUENCE  OF  MODE  CAECUEATIONS 


VELOCITY    FUNCTION 


a  Maximum  Modes   Region   1 

b  Maximum  Modes   Region  2 

e  Maximum  Modes   Region  3 

d  Maximuni  Modes   Region  4 


Figure    19 


84 


APPENDIX  A 
NQRMOl  FINITE  DIFFERENCE  SCHEME 

In  order  to  iterate  equation  (18)  through  the  depth  grid,    a  finite 
difference  scheme  is  used  which  computes  values  of  the  eigenfunction, 
Z(z)  and  its  derivative,    Z'(z)  at  each  grid  point.      This  scheme  was 
originally  developed  by  the  Arthur  D.    Little  Company  [Report  No.    C- 
70673,    31  January  1969]. 

Z(z)  is  expanded  into  a  Taylor  series, 


2!  3? 


where  h  is  the  depth  grid  spacing.      From  equation  (18)  we  have 

.  Z"rz)  =  -  [E-\/ui  Z(z)  =  -  Pu)H(K)  ,  (92) 

where  F(z)  represents  the  square  of  the  vertical  v/ave  number.     Using 
a  forward  difference  for  Z"'(z), 

2'cH)    =    _Z^^(z,U)    -  Z^fz^  (93) 

h 

Substituting  equation  (92),    we  have 

(94) 
Combining  equations   (91),    (92)  and  (94), 

z:  3?L  h  J  *  ^     ' 


h 
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we  have  on  clearing, 

2te.h)ri*l^F(2.h)l    --Zrz)ri-h'F(2)"]    ^  hZ^.  (96) 

Placing  the  above  result  in  a  vector  subscript  notation  we  have  the 
formula  used  to  compute  the  succeeding  eigenfunction, 

We  must  now  compute  the  value  of  the  derivative  of  the  eigenfunction 
Z'(z)  at  the  next  grid  point.     In  order  to  do  this  let  us  express  the 
eigenfunction  at  our  original  grid  point  as  a  backwards  Taylor  series 
of  the  new  grid  point,    z+h, 

V  3?  (Vo) 

Substituting  equations   (92)  and  (93),    we  have 

2.  e 

Rearranging,    we  have, 

ZT.*v.)  -  l[z«-h)(l-f  Rz.h))  -  Z(E)(1  *  I  Fez.))]  ^  (100) 

or  in  vector  notation, 
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APPENDIX  B 
NQRM04  CONNECTION  FORMULAE 

N0RM04  requires  the  evaluation  of  the  phase  at  the  upper  and 
lower  turning  points  in  addition  to  the  phase  effect  of  a  barrier.      For 
the  purposes  of  this  discussion  we  will  align  the  origin  of  the  depth 
axis  with  the  turning  point,    designating  the  positive  direction  towards 
the  "real"  region. 

Consider  the  upper  turning  point  case  (Figure  5).  The  free  sur- 
face boundary  condition  requires  that  the  "imaginary"  region  solution 
disappear  at  the  surface.     From  equation  (59)  this   requires  that 

Z(-H)=K^(-H)[ce^^^-"^  -  De-^^^-^'^]  =0.  (102) 

Thus,    we  can  define  C  and  D  (except  for  a  constant  of  multiplication); 

(103) 
Applying  equation  (71),    we  have 


f°i-v-^>  ::s..H>  •  (104) 


C.--e-^''-"'  D-e^^'""' 


or 


9a  =  Arc  ban  [Txt.  h  (JK^  d^)] 


(105) 
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Now,    consider  the  bottom  boundary  condition  and  lower  turning 
point   (Figure  5).      The  bottom  boundary  condition  (equation  26)  re  - 
quires 


^Kb  =   Ks     Ce!i^"^-  De"^^^-"' 


fb 


Q^^zC-^)  ^    l)^-Sz.i-H)      ' 


(106) 


where  K     is  defined  as  the  imaginary  vertical  wave  number  evaluated 
at  the  water  side  of  the  bottom  interface, 


Ks^=[V(-H)-e]=i,;-_c^^  . 


CC-H) 


:io7) 


Thus,    the  imaginary  solution  can  be  written  by  defining  C  and  D  such 
that 


(108) 


Now,    by  applying  the  turning  point  connection  formulae  (equation  71) 
we  have 


(109) 


or 


0L  =  Arc tan 


D^_ll 


D^  -1 


(110) 
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where  D^  is  defined  as 


(111) 


Now,    consider  the  situation  with  an  upper  and  lower  channel 
separated  by  a  barrier   (Figure  8).     Let  the  z  origin  be  at  the  lower 
connection  point,    z    ,    w^ith  the  positive  direction  downward.     From 
equations   (59),    (69)  and  (71)  we  have  for  a  connection  at  the  upper 
point,    Zi, 


OLi 


De-^ 


Ce^ 


bi       De-"-    -  Ce^ 


(112) 


where  L  is  as  defined  in  equation  (86).      Letting  A-^  equal  the  numer- 
ator and  Bi   the  denominator,    we  have,    except  for  a  constant, 


C=^(ai-k,)e-" 


(113; 


Applying  these  expressions  to  equation  (71),    we  have 


o-i       (ai  +  bi)e^     -^  Cai- bt^e"*- 


or 


b^       (cLt  ^-bi)e'-     -  (ai-bOe""- 


(114) 


[ 


Q^   =  Arcba.1  r(ta.,ei-l)e-.(tan9.-l)e-^" 


(tcLnQi  ♦De'-  -  (toLnQi-l)e 


^ 


(115) 
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APPENDIX  C 
N0RM05  CONNECTION  FORMULAE 

For  N0RM05  the  development  of  the  upper  and  lower  "imaginary' 
region    particular  solutions  by  the  application  of  the  free  surface  and 
bottom  boundary  conditions  is  identical  to  that  for  N0RM04.      The 
differences  in  the  two  programs  lie  in  the  application  of  the  Airy 
phase  connection  formulae  vice  the  asymptotic  connection  fornraulae. 
To  obtain  the  upper  turning  point  phase,    we  apply  equation  (103)  to 
equation  (77) 

"^-JT^^^  ^  (116) 


^  (117) 


or 


0u  -  Arc  tan 


where  D,  is  defined  as 


(118) 


Dx=  Ee^^^^"''^  .  (119) 

Note  that  as  H  approaches  zero  (the  turning  point  approaches  the  sur- 
face) the  value  of  0  in  equation  (119)  approaches  arctan  (1/3)  vice  the 
phase  we  required  for  the  free  surface  boundary  condition,    (0     =  0  ). 
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To  derive  the  lower  turning  point  phase,    we  apply  the  particular 
lower  "imaginary"  region  solution  (equation   108)  to  the  Airy  phase 
connection  formulae   (equation  77), 


t  -  J,  [£  (eb  Ks  -  Co  Kb)  e^-^-'^Lc?bKs*  e^K^)  e"^^'"^]  ^  ( 1 2  i ) 


Dividing,    we  have 

^ED,  ■^   1 


0L  =  Arctafl 


2Da-  1 


(122) 


where  D^  is  defined  as  in  equation  (111).     Note  that  this  equation  also 
does  not  approach  the  bottom  reflected  boundary  condition  (equation 
84)  as  H  approaches  zero  (the  lower  turning  point  approaches  the 
bottom). 

Lauvstad  (1971)  derived  an  expression  for  the  phase  coefficients 
of  a  wave  entering  a  barrier  region  in  terms  of  the  phase  of  the  wave 
leaving  the  opposite  side  of  the  barrier: 

(123) 


75"  L.  _J      • 


By  rearranging,    we  obtain 
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which  agrees  with  equation  (114),    the  N0RM04  result, 


(124) 


0L  ""  Arc  tan 


(ai  -^^1)6^   ^  (ai  -  hi)  e"*" 


(114) 
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APPENDIX  D 


INPUT    DATA 


The  following  format  is  used  for  the  input  data  to  all  three  programs, 
NORMOl,    N0RM04  and  N0RM05. 


VARIABLE 
CARD  ONE 
HED  (20) 
CARD  TWO 

NUMV 


FORMAT         MEANING 


RANGE 


CARD  THREE 


20A4 


110 


NMOD 

no 

lEX 

no 

N 

no 

lOUT 

no 

IDEV* 

12 

H 

F10.3 

CB 

FIO.  3 

DRO 

F10.3 

DRB 

FIO.  3 

HP  LOT 

FIO.  3 

HEADING.    TITLE  CARD 


Number  of  Depth/Sound  Speed    2 
Pairs.     Negative  if  conver- 
sion from  feet  to  meters  de-  ■ 
sired. 


50 


Number  of  Modes  Desired 


Iteration  limit. 


1   -   100 


1   -   14 


Number  of  grid  points  desired    100   -500 

100-1000= 


Type  of  output  desired 
Output  device  for  cardfile 


see  below 


Bottom  Depth 
Bottom  Sound  Speed 
Water  density- 
Bottom  density 
Maximum  depth  to  be  plotted       HPLOT     H 


-N0RM04  and  N0RM05  only. 
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CARD  FOUR 


NFREQ 


110 


Number  of  frequencies 
desired 


CARD  FIVE      (SIX) 


FREQ  (NFREQ)  8F10.3  Frequencies  desired 

CARDS  SIX  TO  5+NUMV     (SEVEN  to  6+NUMV) 


DEP(I) 
CC  (I) 


FIO.  3 
FIO.  3 


Depth 
Sound  Speed 


1   -  20 


DEP(I)> 
DEP(I-i; 


For  NORMOl,    lOUT  =  l   gives  a  punched  deck  output  of  the  sound 
profile  and  mode  parameters  and  profile.     For  N0RM04  and 
N0RM05,    the  lOUT  parameter  gives  the  following  types  of  output; 


lOUT 
0 
1 
2 
3 
4 
5 


PRINTER  GRAPH      PUNCHED  CARDS        DISK  OUTPUT 


YES 

YES 

YES 

no 

no 

no 


no 


YES 


no 


no 


YES 


no 


no 

no 

YES 

YES 

no 

no 


94 


NORMOl  FLOWCHART 


READ  AND 
CHECK  INPUTl 
DATA 


/ 


<7 


LINEARLY  INTERPOLATE 
SOUND    PROFILE 


SOUND  VELOC 
ITY  GRAPH 


PUNCH 

SOUND  VELOC4 

ITY  PROFILE 


FREQUENCY  LOOP  BEGINS    IFREQ^l,  NFREQ 


COMPUTE     VELOCITY 
FUNCTION 


below  cut-off  freq 


£ 


CHECK  FOR 
MAXIMUM  ANDI 
MINIMUM  NUM 
BER  OF  MODES 
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A 


o 


cy 


o 


MODE  LOOP  BEGINS 


M  =  1,    NMOD 


> 


3ZL 


USE  METHOD  OF  HALVES 
TO  GET  UPPER  AND  LOW 
ER  BOUND  OF  EIGENVALUE 
WITH  M-1  AND  M  MODE 
CROSSINGS 


C 


JSZ 


ITERATION  LOOP  BEGINS 


Qim 


) 


sz 


COMPUTE  NEXT  GUESS 
BY      HALVING 


JS 


INCREMENTAL 
CHANGE  <  DESIRED 
ACCURACY? 


no 


yes 


v- 


CALL  ITERAT 


I 


^ 


ITERATIO^ 
COUNTER 


^ 
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m 


SURFACE  MODE 
VALUE  <  ACCURACY 
DESIRED? 


no 


ves 


e^SZ. 


MDDE 


yes^_,,*^  '•^^^  \^  no 

CROSSINGS  <M> 


\7 


P4 
O 

s 

1^ 


PM 
O 

o 
1-^ 


2 


REPLACE   UPPER 
BOUNDS 


8  XZ 


REPLACE   LOWER 
BOUNDS 


<: 


T 


J 


END  OF   ITERATION  LOOP 


) 


10 


CHECK  FOR  DEGENERATE 

SOLUTION  AND  FIXUP 

IF     REQUIRED 


11      4^ 

(     >DDE   HAS    BEEN   FOUNDJ- 


I 


NORMALIZE  MODE   TO  ONE 


^ 


m 
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fi^ 


(14 

8 

p^ 

•-J 

0 

0 

g 

hJ 

0 

g 

P 

o* 

^1 

pj 

Fm 

:^ 


COMPUTE  A  BOTTOM  SOLUTION 


KT 


GRAPH   THE   MODE 


5Z. 


COMPUTE   NORMALIZATION 

INTEGRAL,    GROUP  AND 

PHASE      SPEEDS 


12 


C 


^ 


es 


^> 


PUNCH  MODE 
OUTPUT  DECK 


END  OF  MODE  LOOP 


15 


} 


{ 


END  OF   FREQUENCY   LOOP 


3 


XT' 

(  STOP  J 
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N0RM04  AND  N0RM05  FLOWCHART 


Zi 


READ  AND  CHECK 
INPUT  DATA 


71 


V 


1 


LINEARLY  INTERPOLATE 
THE  SOUND  VELOCITY 


99 


(h 


o 
o 
l-J 

>< 
u 

n 

a- 


i 


COMPUTE  VELOCITY 
FUNCTION  V(Z) 


below  cut-off 


frequency 


P4 
O 

s 


i 


CHECK  FOR  MAXIMUM 
AND  MINIMUM  MODES 


MODE  LOOP 


i 


MODE  =  1,  NMODF 


I 


CALL  SEARCH 

(SEARCH  FOR 
EIGENVALUE) 


i 


CALL  SHAPE4 

C  COMPUTE  THE 
MODE  SHAPE) 


> 


(h(h 


INTEGRATE  MODE,  COMPUTE 
PHASE  AND  GROUP  SPEED, 
AND  DECAY  COEF. 


100 


m 


o 

o 

o 
>*      o 


s 
iii 


PRINT  MODE 
PARAMETERS 


C 


C 


OUTPUT? 


^>^H> 

PUNCH  MODE 
SHAPE  DECK 

1 

>DDE   LOOP  ENDS 


i 


FREQUENCY  LOOP  ENDS 

717 

(  STOP  ) 


} 


) 
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SEARCH  SUBROUTINE  FLOWCHART 


SET    FOR 
SECOND 
DUCT 


no 


(T^ 


SET   SOME 

constants' 


MKE 
FIRST      GUESS 


no 


^ 


rO 


SET  UPPER 
LIMIT 


3  14 


0  0 
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32 


SET  LIMITS 

FOR 
UPPER  WELL 


RELAX 
LOWER 
LINTT 


RE-SET 
LOWER 
LIMIT 


UPPER 
DIFFERENCE   ^  " 


16 


COMPUTE  NEXT  GUESS 
BY  REGITA  FALSI 


17 


CALL  WKB 


I 


ITERATION 
COUNTER 


^^ 


d 


yes 


GO   TO 
TOP  DUCT 


^ 
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16  17 


24 


19 


"A  DIRECT  HIT" 


REPLACE 

UPPER 

LIMITS 


REPLACE 
LOWER 
LIMITS 


FIND  NEXT 
GUESS  BY 
HALVING 


I 


25 


1_1 


THE  EIGEWALUE 

HAS 

BEEN  FOUND 


no^X^  ITERiMIONS 
<  ITMAX? 


ves 


^„ 


26 


1 


27 


NEXT   DUCT 
WILL  BE   CENTER 

RJlCORD   LIMITS 


C 


I 


i 


NEXT  DUCT 
WILL  BE  TOP 

RECORD  LIMITS 


RETURN 


-)     C 


I 


<r  LOWER  DUCT? 
lower  ^s^ 


RECORD 
LIMITS  OF 
CENTER 
DUCT 


RETURN 


3 


^ 
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28 


[ 


I 


CALL  WKB 


yes^y^      IS    THERE 

Ji  LOWER  DUCT? 


I 


NEXT  DUCT 
WILL  BE   TOP 


WILL  LOWER 

DUCT   TAKE  ANOTHE 

MODE 


30     \ 

SET  LOWER 
LIMITS  FOR 
LOWER  DUCT 

V 

NEXT  MODE  WILL 
LOWER  DUCT 

BE 

* 

(      RETURN 

3 

c 


I 


RETURN 


29         V 

NEXT  MODE  WILL 
BE  TOP  DUCT 

1 

(^    RETURN    J 

3 
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WKB  SUBROUTINE  FLOWCHART 


SET  UP  SOME 
INITIAL 
VALUES 


INTEGRATE 

IMAGINARY 

VERTICAL 

WAVE  NUMBER, 


CLOSE  REAL 

VERTICAL 
WAVE  NU^ffiER 
INTEGRATION 
START  IMAGINARY 
WAVE  NUMBER 
INTEGRATION 


INTEGRATE 

VERTICAL 

WAVE  NUMBER 


^-<h 


10 


CLOSE  IMAGINARY 

VERTICAL 

WAVE  NUMBER 

INTEGRATION 


tl 


^ 
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1^ 


12 


£ 


yes 


no 


COMPUTE  UPPER 
PHASE  ANGLE 


15 


I 


yes 


COMPUTE   PHASE 

EFFECT  OF  THE 

BARRIER 


1 


ARE  THE 
TWO  DUCTS 
ONNECTED? 


no 


14 


17 


RE -START   INTEGRATING 
REAL  VERTICAL  WAVE   NUMBER 


1 


ONE  MORE   DUCT 

LOWER   PHASE 
=  'Tf  /  4 


21 


C 


i 


DEPTH  LOOP  ENDS 


25 


I 


) 


BOTTOM 
REFLECTED? 


22 


COMPUTE  BOTTOM  REFLECTED 
PHASE  ANGLE 


1 


COMPUTE  LOWER  PHASE 

ANGLE  FOR  DECAYING 

SOLUTION 


27 


COMPUTE 

DIFFERENCE 

FUNCTION 


C 


<H-L 


RETURN 


) 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c  c 

C       NCRMOl  :  A  NORMAL  MODE  EIGENVALLE  C 

C  AND  EIGENFUN'CTION  SOLVER  C 

c  c 

C       LT  KIRK  E.  EVANS,  USN  C 

c  c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

c 

ir'PLICIT  REAL*8(0) 

c 

C<<<<<<««<<<«<««<<<<<«<«<<<<<<<<<<  COMMON  BLOCKS  >>>> 
C 

CC^MON  /DBLL/  DVZ ( 501 ) , CH, DH2, DCbSG ,DRORB , CKAP, DUMAX 

CC^'MQN  /DdLE2/  DU2  ,  DCB  ,  DC:-1 1  N  ,  Dw  ,  CC  i  ,  DCL  S  t 

COMMON  /D6LE3/  DUZ ( 50  I ) ,OZZ ( bO i ) 

CCNMON  /DtiLE4/  DS  ,  CS5CC  ,  DPV  EL  ,  DC- V  EL 

COMMON  /PARAM/  N,NP1tIEX 

COMMON  /SOUND/  Z(5(.n,C(50J 

CCI^MON  /^lEkO/    H,ED(20) 

CCMMON  /SINbi/  ZM,C8,CMIN,CMAX,PLT^AX,UM,FG 

CCMMGN  /INPUT/  NUMV,NMOD 

CCMMON  /DMAP/  OKI(LOO) 

CCMMON  /DKMAP/  DKM IN  ,DKMAX , DKL , OKU 

CC^f^ON  /DSQUND/  DCC(50i) 

CCMMON  /RHC/  RO,Re 

COMMON  /BGTTO/  DEC  AY, ZB ( 20 )  ,UB ( 20) 

CCNMON  /F^L^/     FQV(20) 

COMMON  /ACC/  CEPS,DIFFI 
C 

C«<<<<<<<«<<««<<«<<<«<<<«<«<<<<<<  DATA  >>>>>»>>»» 
C 

D/^TA  DP2/6.282185306D0/ 

DATA  NPRINT/6/ 
C 

C«<«<««<<««««<<<««<««<<«<<<<  PROGRAM  BEGINS  >>> 
C 


C 
C 


WRITE  (NPRINT,16) 

DC  1  1=1,100 
1  DKU  I)  =-l.CDO 


CALL  INPUTl  (IOUT,NFREC} 

CALL  INTRPL 

DCMIN2=DCMIN=^DCMIN 

CALL    CZPLQT 

IF     (lOUT.NE.l)     GO    TO    2 

CALL    P-^GFIL     (CZZjDCCNPl) 
2    CONTINUE 

NMODP=NMGD 
C 
C  . FRECUENCY    LOOP    EEGINS 

C 

c 

CO    15    IF=1,NFREQ 
FC=FQV( IFJ 
DW  =  0P2^=FQ 

D^2=DW^DW 
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c 

c- 
c 
c 


c 
c 
c 


c 
c 
c 
c 


c 
c 
c 


DKAP=DW/DCNIN 


c 
c 


c 
c 

c 

c 
c 


c 
c 
c 
c 


COMPUTE  THE  V(Z)  FUNCTION 


DC    3 

CCZ=DCC(J) 
CCZ2=CCZ*DCZ 
2    DVZ(J)=DW2=^(  (CCZ-DCMIN)*(DCZ  +  DCMIN))/(DCZ2=*CCMIN2) 


\tiRllE     (NPRII^T,18) 

NNCO=NKCDP 

C/iLL    MAXMIN    (£13) 


FQ 


MCDE  LOOP  BEGINS 


DC  12  M=1,NM0D 

IT=1 

CALL  HALF  (M,£;9) 


ITERATICN  LOOP  EEGINS 


DKN=(DKL+DKU}*0.5D0 
DINC=DABS(CKU-DKL)*C.5DC/DKN 
IF  (DINC-OIFFI)  10,5,5 


5  CALL  ITERAT  CDKM,MSJ 
IT=IT+I 

DEPSIL  =  DUMAX^'^OEPS 

IF  (DA5S(CUZ(NPU  )-DEPSIL) 


11,6,6 


IF  (MS-M)  7,8,8 

DKU=DKI^ 
GC  TO  4 

CKL=DKN 
GC  TO  4 


MODE  HAS  BEEN  FCLND 


S  DKN=(DKL+DKU)*0.5DC 
IC  CALL  FIXUP  (DKN) 

11  DKNAX=OKN 
DK1(H)=DKN 

CALL  DNCRMl  ( DUZ , OUMAX ,NP1 J 

CALL  BOTTOM  (DKN,M) 

CALL  MODPLT  (M,  1) 

WRITE  (NPRINT,17i  IT 

CALL  INTEGR  {DKN,M) 

IF  (lOUT.NE.i)  GO  TO  12 

CALL  OUTPUT  ( DUZ , DKK ,M ,NP1 ) 

12  CONTINUE 


13  DC    14    11=1,100 

14  DKK  II  )=-l.CDO 

15  CC^TINUE 


—  NCDE  LOOP  ENDS 


FREGUENCY  LCCP  ENDS 
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c 

c« 

c 

c 

c 


STOP 
<<<<<<<<<<<<<<««<<<<<<<<«<<<<«<<    FORMAT     ST/iTEMENTS    >> 


it    FORMAT     (•!       K.E.     EVANS    NORMAL    MCDE    E IGE NVA LLE • , 

-'SCLVERt    VERSION     1,     REVISION    6,     CATED    13    ALGUST    73') 

17  FCRMAT     CO  NUMBER    OF     ITERATIONS    =',     15) 

18  FCRMAT     CO  FREQUENCY     :     •,    F6.1,     •    HZ  .  •  ) 

END 
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StBROUTINE  INPUTl  (  IOUT-,NFREQ ) 
IMPLICIT  REAL''^8(D) 

c  * 

C  *  THIS  SUBROUTINE  INPUTS  THE  SCUNC  AND  PLN  DATA, 

C  *  DISPLAYS  IT,  THEN  SETS  SOME  CONSTANTS  FCR 

C  *  LATER  USE. 

C  * 

c 

c 

C<<<<<<<<«<<<<<<<<<<<<<<<<<<«<<<  CONNCN  BLOCKS  »>>»>>>» 
C 

COMMON  /OBLL/  OVZ ( 501 )  ,  DH ,DH2 ,OCeSC ,DRGR 6 , CKAP, DU^AX 

COMMON  /DBLE2/  D W2 , DOB , DCMI N , DW , DC  1 ,DC i SG 

CCNMON  /HEAD/  HEC(20J 

CCHMON  /SINGI/  ZM  ,  CB  ,CKI  N  ,C^'AX  ,  PLTNAX,UM  ,  FC 

CCNKON  /SOUND/  Z(50),C{50) 

CCNMON  /INPUT/  NUNV.NMOD 

COMMON  /PARAM/  N,NPI,IEX 

CCNNON  /DHHH/  DH2S IX , DH2D3 , OHHL F 

CCMMON  /RHC/  R0,k6 

CCMMON  /FREW/  FgV(20) 

CC^^ON  /ACC/  CEPS.DIFFI 
C 

C«<<<<<<<<<<<<<<<<«<<<<<<<«<<<<«<<<    DATA    >>>>>>>>»>>>» 
C 

CATA    NREAD,NPRINT/5,6/ 
C 

C<<<<<<<<<«<<<«««<<<«<««<<«<<<<    PROGRAM    EEGINS    >>>» 
C 
C 

CCNVRT=l.OEO 
C 
CI ■ READ    IN    TFE    RUN    PARAMETERS 


C 


RE/SD    {NREAC,4)    HED 
WRITE     (NPRINT,5)     HED 


R£^D    (NREAD,6)    NUMV,NJMOD,  I  EX,  N,  I  CUT 

N=IABS(N) 

IF    (N.GT.500)     N=50C 

IF    (N.LT.5G)    N=LOO 

NNCD=IAES(NMOD} 

IF  (NMGD.GT.IOO)  NM0D=10G 

IF  (NMQC.EQ.OI  NM0C=10 

IF  (N.EC.OJ  N=500 

IEX=IABS( lEX) 

IF  (IEX.GT.I4)  IEX=14 

IF  (  lEX.EQ.O)  IEX  =  7 

NFl=N+l 

RE/^D  (NREAD,7)  ZM  ,  CB,  RO  ,RB,  PLTMAX 

2M=ABS(ZM) 

PLTMAX=AdS(PLTMAX) 

RE/iD    (NREAD,6)     NFREC 

NFFEQ=IABS(NFREQ) 

IF  (NFREC.GT.20}  NFPEQ=20 

RE/iD  (NREA0,7)  {  FG  V  ( I  )  ,  I  =  i  ,  NFRE  C  ) 

F£=F&V( U 
C 

C— ' r ^-CONVERT  FPCM  FEET  TC  METERS 

C  IF  REQUIRED 


C 


IF  (NUMV)  1,2,2 

eCNVPT=0.3048 

NUMV=-NUMV 

ZM=ZM*COMVRT 

Ce=CB-CCNVRT 

PLTMAX=PLTMAX*CONVRT 

y^RITE  (NPRINT,8) 

V^RITE  {NPkI(MT,9)  (FGV(  I)  ,I  =  1,NFR-EC) 
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VnRITe  {^PRI^iT,lO)  nnociex 

ksPITE  (NPRINT.Il)  ZN'.PLTMAX 

WPITE  (NPRINT,L2)  POtRB 

WRITE  (NPRINT,L3i  CB,N 
C 

c 

C READ    IN    /»KC    DISPLAY    SCUND    DATA 


WRITE     (NPRINT,14)     HED 

DC    3    I=l,NUMV 

READ    (NRcAC,7)    Z( I ),CI 

CI=CI^CCNVRT 

Z( I)=Z( IJ^CGNVRT 

C(I)=CI 

3  WRITE    (NPRINT,15)     Z(I) ,CI 
C 

c 

Q COiMPUTE    SCNE    CONSTANTS 

C 

Dh  =  CBLE(ZMi/DFLOAT(i\) 

D2^'=D6LE(■ZM) 

CRCRb=DeLE(RO/RB) 

Ch2  =  0ri-'CH 

CCe=DBLE(Ce) 

CCBSQ=CCe^CCB 

DI-2SIX=DH2^0.1666  66b66  6666  667D0 

D(-2D3=Dh  2*0.3  3333  33  33  3333  33300 

DhhLF  =  Dh-«0.5D0 

DEPS=l.COr^^M-IEX) 

CIFFI  =  DEPS-'<DEPS 

CIFFI  =  CNA)^i(DIFFI,1.0D-L6) 
C 

RETURN 
C 
C 

C«<<<<<<«<<<<<<«<«<<<<<<<«<<<«<<<    FORMAT    STATEMENTS    » 
C 

c 
c 

4  FCRf^AT     (20A4) 

5  FCRNAT  {• !• »20A4J 
e    FCRMAT  (5110) 

?  FCRMAT  iSFlO.3) 

8  FCRMAT  COS  T34,  'RUN  PARAMETERS') 

9  FCRMAT  CO' ,T30, 'FRECUcNCIES  =',  (T42,F6.l,«  HZ'»/)) 
IC  FCRMAT  (•0'iT24,   'MODES  RECUESTEC  :'t  14,  //, 

1  T24,  •  ITERATION  LIMIT  =  LO-'t-'C  -',12,')') 

11  FCRMAT  (•C',T27,  'BCTTOM  DEPTH  :',Fe.l,'  METERS',//, 

1  T22,  'MAX  DEPTH  PLOTTED  :'tF6.1,'  METERS') 

12  FCRMAT  ('a',T2o,  ' CCEAN  DENSITY  :',  F7.4, 

*  •  GM/CM^^^3',  //, 

1  T25,  'BOTTOM  DENSITY  :',F7.4,'  GM/CM**3«) 

13  FORMAT  CO'tTie,  'BCTTCM  SOLND  VELCCITY  :', 

*  F8.1,:'  METERS/SEC',//, 

1  T17,  'NUMBER  OF  GRIC  POINTS  :',  15) 

1^  FCRMAT  CI',  20A4,  //,  T33 ,  'VELCCITY  PRCF  I  LE ' , // » T22 , 

1  'DEPTH,  METERS' ,T44, 'SOUND  VELOCI TY  ,M/S£C '  ,  / ) 

15  FCRMAT  CO',  T25,  F6.U  T57,  F6.1  ) 
END 
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SLEROUTINE  IKTRPL 
IMPLICIT  REAL*8(D) 
C 

c  ^ 

C  *  THIS  SUBROUTINE  INTERPOLATES  THE  DEPTH  /SND  SOUND 

C  *  VELOCITY  DATA  GIVEN,  FOR  THE  GRID  VALUES. 

C  *  IN  ADDITION,  IT  CCi.vPUTES  SOi^E  CCKSTANTS  FOR 

C  *  LATER  USE,  AND  EXTRAPOLATES  THE  SOUND  VELOCITY 

C  *  TO  THE  BOTTOM  IF  RECUIRED. 

C  * 

CCFMCN  /DBLl/  DVZ ( 501) , CH , DH2 , DCESC  ,DRORE , CKAP2, DUMAX 

CCMN.ON  /DdLEZ/  D^.2  ,  CCB  ,  DCMI  K  ,DV.  ,CC  I  ,DC  ISC 

CCNKON  /DDLE3/  DUZ  (  501  )  , DZZ ( 501 ) 

CCMMQN  /PARAM/  N,NP1,IEX 

CCf^K:ON  /SINGI/  ZM,CB,CMIN,CMAX,PLTNAX,UM,FG 

CCNNON  /USCUNC/  DCC(p01J 

CCMMON  /SOUND/  Z(5C),C(50) 

CCNMON  /INPUT/  NUMV,NMOC 


C 

c 


CNAX=0.i 
CMN=CB 


C 

C      ***  EXTRAPOLATE  THE  PROFILE  TO  THE  BQTTO^^   *** 
C 

IF  (Z(NUMV)-^lM)  1,2,2 
C 

1  NLNV=NUf^V+l 
C 

2  C(NUMV)=C(NUMV-1}H-(ZW-Z(NUMV-1)  )*C.  17E' I 
Z(KUMV) =ZM 

C 

3  NUNVPI=NUMV+1 
NL^VD2  =  KUr<V/2 

c 
c 

C«<<«<««<«<«<<«<<<«<<  SWAPPING  AND  TESTING  LOOP 
C 

CC  4  I=l,NUMVD-2 
NLF=NUMVP1^I 
CL=C( I) 
CU=C(NUP) 
ZL=Z(I) 
ZL=Z(NUP) 

CNAX=AMAXl(CMAX,CL,CU) 
CNIN  =  A.MIM{eMIN,CL,CU) 
C 

C(I)=CU 
C(NUP)=CL 

z(n=ZM-zu 

Z(NUP)=2M-ZL 
A    CONTINUE 
C 

C<<<<«<<<«<<<<«<«<<<<<<<<  SWAPPING  AND  TESTING  LOOP  ENDS 
C 

IF  <NUMVD2-^2.EQ.NUMV)  GC  TO  5 

CNVD21  =  C(NJMVD2-^l  ) 

CyAX=AMAXl (CMAX,CNVC21 ) 

CNIi\=A:-lIia{CMINTCNVDZl) 

Z(NUiMVD2+i  )  =  Z.M-Z(AUMVD2+1) 
C 
C  ***      SET    UP    FQR    T>iE    FINE    GRID    CONFUTATION       ^** 


C 


DCyiN=DBLE(CMIN) 

CCMN2=LCf^.IN^-DCMIT^ 

DC1=0BLE(C{1) ) 

CCiSQ=CCl^-CCl 

C2^.=0BLE(ZM) 

1  =  1 

CCI=DC1 
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DCIP1=CELE(C(2) ) 
C2I=DBLE(Z(U  J 
C2IPl=CeLE(Z(2) ) 
D2DIFF=DZIF1-DZI 
CFLN=DFLGAT(N ) 
CZF=C.OCO 
C 

c 

C«<««<«««<<«««<<«<<  FINE  GRID  LOOP  BEGINS 
C 

CC  8  J=l,NPi 

C22( JJ=CZP 

IF  (OZP-DZIPL)  7,7,6 
C 

C  ***   RESET  THE  I  PARAMETERS   *'** 

C 

6  1=1+1 
D2I=DZIFI 
CCI=DCIPI 

CCIPI  =  CELE(C(  I+U  ) 

CZIPI=DELE(Z(H-l)  ) 

C2CiFF=DZIPl-DZI 
C 

C  «**   COMPUTE  THE  VALUES  FOR  C(Z)   *** 

C 

7  DC2=(DCI^{DZIP1-DZP)+DCIP1*(CZP-CZI ))/DZDIFF 
CCC(J)=DCZ 

c 

C         ***   COMPUTE  THE  NEXT  DEPTH   *** 
C 

C2P  =  DZM='^DFLnAT(J  iVDFLN 

8  CONTINUE 
C 

C<<<<<<««<<<<<<««<<<<<<<<  FINE  GRID  LOOP  ENDS 

C 

C 

RETURN 

ENC 
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SLBHOUTINE  MAXMIN  (*) 
IMPLICIT  REAL*8(D) 
C 

c  * 

C  *  THIS  SUBROUTINE  CHECKS  THE  NAXIMUM  ANC  MINIMUM 

C  *  VALUES  OF  THE  HCRIZCiMTAL  WAVE  NUMBER;  ThEN 

C  *  DETERMINES  THE  NUMBER  OF  MODES   REPRESENTED 

C  *  BY  EACH.   IF  NC  MODES  ARE  PRESENT  FOR  THE 

C  *  MINIMUM  WAVE  NUMBER  THE  FRECLENCY  IS  EELCW 

C  *  CUT-CFF.   IF  THE  NLMdER  OF  MGCES  REQUESTED  IS 

C  *  ABOVE  THAT  REPRESENTED  BY  TFE  MINIMUM  r^AVE 

C  *  NUMBER,  THE  REQUEST  IS  MODIFIED  ACCORDINGLY. 

C  * 

c 

C 

C<<<<<<<«<<<<<<<<<«<«<«<<<<«<«<<<  COMMON  BLCCKS  >>>»> 

C 

C 

CCMMON  /DBLl/  DVZ ( 501 ) ,DH ,0H2 , DCBSC ,DRORB » CKAP , DUMAX 

CCMMON  /DBLE2/  DW2  t  CCB  ,  CCM  I  N  t  OW  ,  CC  1,DC  LS«Q 

COMMON  /DBLE3/  DU Z (501 ), DZZ ( 50 1  ) 

CCMMON  /PARAM/  N»NP1,IEX 

CCMMON  /HEAD/  HED(20) 

CCMMON  /SINGl/  ZM , C3 ,CM IN ,CMAX » PLTNAX ,UM , FG 

CCMMON  /CMAP/  DKL(IOO) 

CCMMON  /DKMAP/  DKM  IN ,OKMAX » CKL , CKU 

CCMMON  /INPUT/  NUMV,NMGD 

CCMMON  /SOUND/  Z(50),C{50) 
C 

C<<<<<<<<<<<<<<<<<<«<<<<<<<<«<<<«<<<  DATA  >>>>>>>>>>>>>» 
C 

DATA  NREAD,NPRINT/5,6/ 

DATA  MAXI,MINI/«MAXI« ,  'MINI'/ 
C 

C<<<<<<<<<<<<<<«<<«<<<<<«<«««<<<<  FROGRAN  BEGINS  »>» 
C 
C 

C      **«   COMPUTE  THE  MAXIMUM  AND  MINIMUM  DK'S   *** 
C 

CKMN  =  OW/DCB 

DKMAX  =  0^^/DCMIN 
C 

C      ***   CHECK  OUT  THE   OKMIN  FOR  MAX  NO.  OF  MCCES   *** 
C 

CALL  ITERAT  (DKMIN,MrS) 

IF  (MS)  1,4,1 

1  MCCMAX=MS 
DK1(MS)*DKMIN 

C 

MCCMIN=0 

MS  =  0 
C 

C      ***   PRINT  THE  RESULTS   *** 
C 

WRITE  (NPRINT,5)  HED 

WRITE  (NPRINT,6) 

WRITE  (NPRINT,7)  M  AXI  ,  CKM  AX  ,  MOD  MN 

WRITE  (NPRINT,7)  f^  INI  ,  OKMI  N  ,  MOD  MAX 
C 

C      ***   CHECK  THE   NUMBER  CF  MCDES  PECUESTEC   *'^* 
C 

IF  (MCCNAX-NMOD)  2,3,3 

2  NNCD=MOOMAX 
WRITE  (NPRINT,3) 
WRITE  (N PR  I  NT, 9) 
WRITE  (NPRINT,10) 
WRITE  (NPRINT,9) 
WRITE  (NPRINT,8) 
WRITE  (NPRINT,il)  NMOD 
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2  RETURN 

*=*»?   6ELGW  CUT-OFF  FREQUENCY   **'5= 


UPITE 
WRITE 
WRITE 
WRITE 
WRITE 
WRITE 
RETUkN 


(NPRINT,8) 
(NPRINT, 9) 
(NPRINT, LO) 
(NPRINT,9) 
(NPRINT,8) 
{NFRINT,12) 
1 


C 

c 

C<<<<<<<<<<«<<<<««<<<<<<<<<«<«<<<<    FORMAT    ST/iTEMENTS    » 
C 

c 
c 

c 


8 
c 

10 
11 
12 


5  FCR 

6  FCR 
1 

7  FCR 
1 

FCR 
FCR 

FCR 
FCR 
FOP 


NAT 
MAT 

MAT 


(•IS 


CN  HORIZONTAL* 
,  F12.8,  3X, 


20A4i 
(•0',  ///, 

•  WAVE 
CO',  T22t 
•MODES  :' 
(T33,13( '* 
(T33, 1H-,11X,1H*) 
(T33t'=?=   CAUTION   *') 

(•O't  T26,  'MODES  RECUESTEC  RESET  TO  : 
CO',  /,  T25,  'FREQUENCY  RECUESTEC  IS 
CUT-OFF',  //,  T20,  'JOB  TERMINATED 


T25,  'BOUNDS 

NUMBER'  ) 
A4,  'NUM  K  : 

}) 


ENC 


NAT 

MAT 
MAT 
MAT 
MAT 

•  BELOW 

•NEW  FREQUENCY  SELECTED.') 


14) 


OR 
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SLEROUTINE  HALF  (M,*) 
I^FLICIT  REAL-8(C) 
C 

Q  *:»:k3|t^:>;j?3h*«j;i5k**3{c**«3|.  ****♦*******:«!**************  ^:  ******* 

C         ^ 

C  *  THIS  SUBRCUTINE  SEARCHES  FOR  AN  UPPER  VALUE 

C  *  OF  CK  WITH  M-l  KODE  CROSSINGS  A,\D  A  LG^nER 

C  *  VALUE  CF  DK  WITH  M  f^GUE  CROSSINGS.   AT  THE 

C  *  SAME  TIME  IT  FILLS  IN  A  MAP  CF   DK  VS. 

C  *  MODE  CROSSINGS  FOR  LATER  USE. 

C  * 

Q  5{=***********'':**^***  *>'****  A  **********■*-*****  ************ 

C 

C«<<<<<<<<<<<«<<<<<<<<<«<<<<<<CCMNCN  BLOCKS  >>>>>>>>>>>>> 

CCNMCN  /DBLE2/  DW2 » CCB t COM IN,DW  ,CC  I  ,DC LSC 

CCNFON  /PARAM/  N,NPI,IEX 

CCMNON  /DMAP/  DKI(ICO) 

CCMMGN  /DKMAP/  DKM IN, DKMAX , CKL , DKU 

CCM^ON  /ACC/  DEPStCIFFI 
C 

C<<<<<<<<<<<<<<<<««<<<<<<<<<«<<  PROGRAM  BEGINS  >»»»>>> 
C 

LCW=-l 

LUf=-l 
C 
C      *'?'*  CHECK  THE  MAP  FOR   ALREADY  COMPUTED  GLESES   *** 


C 

c 


c 
c 
c 


c 


c 


DC  2  I=Mt 100 

IF  (DKl(I)  )  1,1,3 

1  CONTINUE 

2  CONTINUE 


DKL=DKMIN 
GC  TO  4 


3  DKL=DK1(  1) 

4  CKU=DKMAX 
C 

C 

C      ***   FIND  £  CHECK  THE  NEW  MID  DK  VALUE  *** 

C 

5  Ci<N  =  (DKL+DKUJ*0.5DO 
C 

C      ***   CHECK  THE  DIFFERENCE   ** 


IF  {(  CKU-DKN)/DKN-DIrFn  13,6,6 
e    CALL  ITERAT  (DKN,MS) 
IF  (MS-N)  7,8,9 


C 

C      ***   IF  WE  GET  TWO  END  POINTS  WE  ARE  FINISHED   *** 
C  ***   IF  NCT,  CONTINUE  HALVING   *=«■* 


CKU=OKN 

IF  (LUP.GT.O)  RETURN 

LCW  =  1 

GC  TO  5 


8  CKL=DKN 

iF  (LOW.GT.O)  RETURN 

LUF=1 

GC  TO  5 
C 
C      ***   FILL  IN  THE  MAP  IF  NEEDED   *** 


9  IF  (MS- ICC)  10, 10,£ 
10  IF  (CKN-DKKMS)  )  12,12,11 
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c 
c 
c 


c 

c 


11  CK1(MS)=DKN 

12  CKL=DKN 
GC  TO  5 

^^^       IF  WE  GET  HERE  WE  HAVE  COME  TC  WITHIN  *»*« 
***   THE  LIMITS  OF  DESIRED  ACCURACY   KM  DK  *** 

13  CKL=DKK 
CKL=OKN 
RETURN  I 


END 
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SLBROUTINE    HERAT     (DKtNS) 
I.VPL1CIT    REAL-<3(D) 
C 

c 

C  * 

c  * 

C  *  THIS  SUBROUTINE  ITERATES  THE  FINITE  DIFFERENCE 

C  *  EQUATION  FROM  THE  BOTTOM  TO  THE  SURFACE.   THE 

C  *  COMPUTED  VARIABLES  ARE: 

C  * 

C  *        DiJZ  (I)  :  VALUE  OF  THE  NCCE  E I GENFL  ^CTI  ON 

C  *        ^^S  :  NUMBER  OF  MODE  CROSSINGS 

C  *        DUMAX  =  ABSOLUTE  MAXLMUN  VALUE  CF  EIGENFUNCTIO 

c 
c 

C<<<<<<<<<<<<<<<<<<<«<<<<<<<<«<«<<<<    COMiMON    ELCCKS    >>>>» 

C 

C 

CCNMCN  /DELI/  DVZ ( 501 ) t CH, DH2 , DCBSG ,DRORB , DKAP , DUMAX 
CCKMON  /DBLE2/  0W2  ,  CC  B  ,  COM  1 1\  »  DW  ,CC  L  »DC1S  G 
CCMMON  /DBLE3/  DUZ ( 50  I ) ,OZZ ( 50i ) 
CCNFCN  /PARAM/  N,NF1,IEX 

CCF.MON    /SINGl/    ZM,CB,C>aN,Cf^AX,PLTNAX,UM,FG 
CCMf/ON    /CKf^AP/    DKMIi\|,CKMAXt  CKL,CKL 
CC/^MON    /OHHH/    DH2S  IX  t  DHzD3  »  CHHL  F 
C 

c 

C«<<<<<<<<«<<<<<<<<<<<<«<<<<<<««<<  DATA  >>>>>>>>>>>>>>> 
C 

DATA  CUPPER/C  .1D6C/ 

DATA    NPPIiNT/6/ 
C 

C<<<<<<<<«<<<<<«<<<<<<<<<<<<«<«<<<<  PROGRAM  EEGINS  >>>» 
C 

c 

C . DEFINE  SCNE  CONSTANTS 

C 

DSTART=i.CC-10 

CK2=DK-CK 

CE=(DKAP-OK)*{DKAP+CK) 
C 

C- r — ' SET  UP  RUN  PARAMETERS 

C 

1  CLVAX=OSTART 

CN>lRK=i.OD0 

NS  =  0 
C 

C ' INITIAL  VALUES 

C 

CUJ=DSTART 

CLFJ=DkCR6«DSCRT  (  (  CK-CKMIN  )*  (  DK■^DK^'IN)  )*DSTART 

CFJ=D£-CVZ(1) 

CLZ(1)=CUJ 

c 

C  — ' ^ -— ^ — DEPTH  LOOP  BEGINS 

C 

c 

DC  8  J=2,NPI 
C 

C ITERATE  THE  FUNCTION 

C 

CFJPL=LE-DVZ{J) 

CENCM=l.0DC+DH2SIX*CFJPI 

DLJP1  =  (  {I.0D0-DHzD3^CFj  J>^:=CUJ  +  Dh=?CUPJ)/DENCN 

D^^(  l.GGC-DH2C3^DFJPl  )^:=CUPJ 

CE=(i:FJ-fDFJPl-DH2S  IX^  C  F  J-<OF  JP  U  =*DU  J 

CLFJP1=(CA-D5^'DHHLF)/DEN0M 
C 
C . ^-CHECK  FOR  MAX  DU(J) 

C 
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CAUJP1  =  C/^BS(DUJPI) 
IF  (OAUJPl-OUMAX)  4,^,2 
2  IF  (DAUJFl-LUPPER)  3,3,9 
2  CLNAX=DA0JP1 
GC  TO  7 
C 
C CHECK  ZERC  CROSSINGS 

C 

4  IF  (DUJF1*DMARK)  5,6,7 

5  ^S  =  F.S+I 
C^ARK=DSIGN{0MARK,CUJP1} 
GC  TO  7 

C 

e    IF  (DUJPl)  7,5,7 

c 

Q -^ ^ SET  POR  (vjEXT  ITERATICN 

C 

7  CFJ=DFJP1 

CLJ=DUJF1 

DLFJ=DUPJP1 

DLZ{ J)=DUJ 
e  CCNTINUE 

c 
c 

Q ^ ._-. CEPTh  LOCP  ENDS 


RETURN 


C 

c 

C OVERFLOW  PROTECT  I  CN 

C  REDUCE  SCALE  OF  U(Z) 

9  IF  (DSTART .LT.L.OD-60)  GO  TC  10 


DSTART=DSTART*i.OD-IO 
GC  TO  I 


IC  V^PITE  (NPRINT,  ill  CK 
^'S  =  C 
RETURN 


C 

C<<<<<<<<<<«««««<<<<<<<<<<««<<<<    FORMAT    ST/iTEMENTS    » 

C 

11  FCPVAT  (•!  SUBROUTINE  ITERAT  OVERFLOW  PRCTECTICN*, 

*  •  EXCEEDED.  ',  //, 

1  '   NS  SET  TO  ZERO  ^^ND  RETURNEC, 

*  •  TO  CALLING  PRCGRAM.',  //, 

2  •    DK  =  '  ,  G2G.12  ) 
ENC 
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SLEROUTINE  FIXUP  (CKN) 
INFLICIT  REAL*3(D) 
C 

c 

c  * 

C  "*  THIS  SUBSOUTIWe  FORCES  A  CECAYING  E)(PCMdNTIAL 

C  *  ON  THE  EIGENFUNCTIGN  A3CVE  TFE  UPPER  TURNING 

C  *  POINT. 

C  * 

c 
c 

C<<<<<<<<<<<<<<<<<«<<<<<<<<<^«<<«<<<  CCMMON  BLCCKS  >>>>» 
C 

CCNf/GN  /PARAM/  NtNFI,I£X 

CCFMGN  /D3LL/  DVZ ( 50  I )  , CH ,DH2 , DCBSC  ,DR0R3  ,DKAP, DUMAX 

CCM^QN  /CaLE3/  DU 2 (501 ),  CZZ (50 i  ) 
-CCMMON  /DHHH/  DH2S  I  X , Dri203 , DHHLF 
C 

C<<<<<<<<<<<<<<<«<«<<<<«<<<<<<«<<<<  PROGRAM  EEGINS  >>>» 
C 
C 

DE=(OKAP-OKN)*(DKAP+DKN)   - 

NF2=NPl+l 

CAPG=DE-LVZ(NP1) 

IF  (DARG.GT.O.ODG)  RETURN 

IFASS=^l 

IF  (DABS(DUZ(NP1 )  )  .LT.DUMAX}  IPASS  =  i 
C 

C- FIND  THE  UPPER  TURNING  POINT 

C 
C  . 

DC  1  1=2, NPl 

INDEX=NP2-I 

CKZ2=0E-DVZ( INDEX) 

IF  (DKZ2)  1,2,2 

1  CCNTINUE 
C 

RETURN 
C 

C      .— ^ DECAY  THE  SCLUTIGN 

C 

2  JTFU=INCEX 
JTPUP1==JTPU+1 
ISTCP=NP1-JTPU 
CUJ=C.OCO 
CLFJ=I.CD-5 
CFJ  =  D£-DV»Z{NP1) 

C 

CC  3  J=1,IST0P 

INDEX=NP2-J 

CUZ(INC£X)=DUJ 

CFJP1  =  D£-DVZ(  INDEX-l) 

CENOM=L  .0D0+DH2SIX'-^DFJPl 

DUJP1=(  (i.0D0-DH2C3^CFJ)*CUJ  +  DH=*CUPJ)/DENCM 

C/=(  1.0CC-DH2D3'*0FJP1  )=^-DUPJ 

CE=(CFj+uFJPI-DH2SIX=^CFJ=pDFJPL)*DLJ 

DUFJPl  =  (DA-Cb'i*DHHLF)/CENOM 

CFJ=DFJP1 

CLJ=OUJFI 

3  DLPJ=DUPJP1 


C 
C 

c 
c 


CCCEF  =  DUZ(  JTPD/DUJPI 
CC    4    I=JTPUPl,N 

A  Duz(  i)  =  cuz(niocaEF 

IF     (IPASS.GT.O)     RETURN 
DUMAX=G.ODO 

CC    6    I=1,JTPU 
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c 
c 
c 


C;iCUZI  =  CABS(OUZ(  I  )  ) 

IF  (OACUZI-DUMAX)  6,6,5 

5  CLNAX=UADUZI 

6  CCMINUE 


RETURN 
END 
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c 
c 
c 
c 
c 
c 
c 

c 

c 

c 
c 


c 
c 
c 


SteROUTINE  CNCRMl  (CLZ , CUMAX,NP 1 ) 
INFLICIT  REAL'^8(DJ 

* 

*  THIS  SLBRCUTINE  NGRNALIZES  THE   DOUBLE 

*  PRECISICN  VECTCR  Lbl     (NPl),  TC  THE  VALLE 

*  DUMAX. 

DINENSICN  DUZ(NPl) 


C^EST  =  DlJ^^AX-0.1D-60 


CC  1  1=1, NFl 

DLZI  =  DUZ{  I  ) 

IF  (DAaS(CUZI) .LT.CTEST)  DUZI=O.0D0 
1  DLZ(I)=DUZI/DU^1AX 


RETURN 
END 
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SLEROUTINE  BOTTOM  (DK,I^) 
INFLICIT  REAL^8(D) 

c  * 

C      *        This  SUBRCL'TINE  COMPUTES  A    60TTGN  SCLND 
C      *    PRCFILE  FOR  PLOTTING  PURPOSES. 
C      =^ 

c 

C<<<<<<<<<<<<<<<<<<<«<<<«<<<<<<<«<<<    COMMON    BLCCKS    >>>>» 
C 

CC^^'CN    /DBLl/    CVZ(501  »  ,CH,CI-2TDCeSCiDR0RBtCKAP,DUMAX 

CCMMON    /DBLE2/    0'/^2  ,  DCB  .  DCM  I  N  , D W  tCCi  ,0C1  SG 

CCf^MON    /DELE3/    DUZ  (  501  )  ,  DZZ  (  501  ) 

CC^.XON    /SINGl/    Z.M,CB,CMIN,CNAX,PLTNAX,UM,FC 

CCKMON  /BOTTO/  DECAY, ZB ( 20  J  ,UB ( 2C ) 

C 

C«<<<<<<<<«<<<<««<<<<<<<<<<<<<«<<<  FRQGRA^'  EEGINS  >>>>> 

C 

Q COMPUTE  CEPTh  VALUES 

C 
C 

HE=(ZM-PLTyAX J*0.5E-i 

DhB=DBL£(F6) 
C 

C J COMPUTE  SOME  CONSTANTS 

C 

CGAM=DSQRT(  (  DK-DW/DCE  ) '!' (  DK  +  DW/DCB  )  ) 

DARG  =  OGAM--CHB 

IF  (DtXFU+OARG)  3,3,1 

1  CEX=CEXP(DARGi 
CV  =  DUZ(  1)-CKCR6 

C 

C : COMPUTE  THE  PROFILE 

C 

c 

CC  2  1=1,20 

IF  (DV.LT.l. 000-30)  GC  TO  4 

DV-DV'^DEX 

2  UE(I)=SNGL{DV  ) 
C 

RETURN 
C 
C- . . I^UL L  PR 0 F  I L E 

C 

3  1=1 
C 

4  CC  5  J=I,2C 

5  UE(J)=0.0 
C 

RETURN 
ENC 
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SLERCUTINE  INTEGR  (CK,NCDE) 
IMPLICIT  REAL^^8(D) 

c     * 

C  *  THIS  SUBROUTINE  COf^PLTES  THE  NC  RI^AL  I  Z/^T  I CN 

C  *  INTEGRAL,  AND  THE  SCLND  VELOCITY  NORMAL- 

C  *  IZATICN  INTEGRAL.   THE.M  THE  GRCUP  AND  PHASE 

C  *  VELOCITIES  ARE  CCMPUTED. 

C  * 

C 
C 

C«<<<<<<<<<<<<<<<<«<<<<<<<<<<<<<<<<<<  COMMON  BLOCKS  >>>>» 

c 
c 

CCN^.ON  /PARAM/  N,NPI,IEX 

CCNMCN  /SINGI/  ZM  ,  CB  ,CM I N , CNAX , PLTNAX ,UM , FC 

CCMMON  /DEL!/  OV^Z  (  50  I ),  LH  ,  DH2  ,  DCdSG  tORORE  ,  CKAP  ,  DUMAX 

CCNKON  /DBLEZ/  DW2  ,  CCB  ,  LCMl  N  ,  DW  ,  CC  1,  DCl  S<* 

CCNMON  /D6LE3/  DUZ ( 5C1 ) , DZ Z (50 1 ) 

CCK>!ON  /DELE4/  DS  t  CSSCC  ,  DP  V  EL  ,  DGVEL 

CCf^CN  /PHC/  RO,RB 

CC^'^'ON  /OSCUND/  DCC(501) 

CCNNCN  /BCTTO/  DEC  AY, Zb ( 20 ) t U8 ( 2C) 

cc^^'ON  /cknap/  dk^^in,ckmax,ckl,cku 

ECLIVALENCE  (DUZ  1 , CLZ ( L ) ) 

c 

C<<<<<<<«<<<«<<<<<«<<<<<<<<<<<<<<<<<<  DATA  >>>>>>>>>>>>>> 
C 

CnA    NREAC,NPRINT/5,6/ 
C 

C<<<<<<<<<<«<<<«<«<<<<<<<<<<<<«<<<<  PRCGRA^<  BEGINS  >>>» 
C 

c 

C  ***         SET    UP    THE     END    POINTS    FOR    THE    TRAPAZCIDAL    RULE 

C 

DLZlSC^CUZd  ) 

IF     (DABS(DLZLSQ) .LT. l.CC-30)    DUZlSC=O.ODG 

CLZlSQ=CUZiS(.«DUZlSC 

CS  =  C.5DC-(DUZ1S>^) 

DCISQ=CCC( 1)*DCC( I ) 

CLZ1SG=CLZ1SQ/(CC1SC) 

CSSCC=G.5D0«(DUZ1SG) 
C 
Q INTEGRATING  LOOP  BEGINS 

C 
C 

CC  1  1=2, N 

CLZI2=CLZ( I) 

IF  (DAES(DUZI2) .LT.i.CD-30)  GO  TC  1 
-.  CLZI2  =  CUZI2--i'CUZI2 

CS=DS+0LZI2 

CCCI=DCC(I) 

ClZI2=CyZI2/(CCCI*CCCII 

DSSCC=DSSCC+DUZI2 
I  CONTINUE 
C 

c 

C INTEGRATING  LOOP  ENDS 

C 

c 


I 


C      ***    ADD  BOTTOM  EXPONENTIAL  DECAY    *** 
C 

CA2  =  (DK-CKMN)*(CK  +  CKMN) 

CE  =  CbLE(R6)=^DUZl-C.50C/DSQRT(DA2) 
CS=CELE(RCi*OH*DS+Ce 
CSSCC  =  DSSCC=^DBL£(  RC}*DH  +  DB/CCBS  C 
CECAY  =  CkN*CB/(CCB»=DK*DSJ 
C 

c 

C      **♦   CONPUTE  THE  PHASE  AND  GROUP  VELOCITIES   *** 
C 


125 


c 
c 


CFVEL=Cl\/DK 
OGVEL=CS/(CFVEL^DSSCC) 

k^PITE     (NPRINT,2)     ^'GCE,DPVEL 
WRITE     (NPRINT,3)     NCDE,DGVEL 


RETURN 


C 

c 

C<<<<<<<<<<«<<«<<<<<<<<<<<<<<<<<«<<<    FORMAT    ST/iTEMENTS    » 

C 

C 

2    FCP^^AT     COPHASE    VELCCITY    FOR    MOCES     14,     •     IS     •, 

i  G15.7t     •     f'/SEC) 

2    FCRr'iAT     COGRGUP    VELCCITY    FOR    MOCES     14,     •     IS     •, 

1  G15.7,     •     N/SEC) 

ENC 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c« 

c 


SLEROUTINE  MODPLT  (M,JUMP) 
INFLiClT  RbAL^8(D) 


* 

THIS  SlibRQUTIKE  PLCTS  THE  MCCE  SHAPE 

CN  THE  PRINTER  USING  THE  N.P.S.  ROUTINE 

UTPLOT  . 


IF    JUMP=1,     ALL    POINTS    IN    THE    MODE 
FUNCTION     WILL    EE    PLCTTEu.       IF    JUNF 
IS    NOT    EQUAL    TO    I,     EVERY    NPl/IOO    TH 
POINT    WILL    BE    PLOTTED. 


<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<«<<<    CCMi^iON    BLCCKS    >>>>>» 

CCNMCN  /DBLl/    DVZCSOU ,CH,DH2,DCESC,CR0RB  ,CKAP2,DUNAX 

CCNMGN  /DbLE2/    DV^2  ,  CCB,  CCM  IN,  DVv  ,  CCl  tOC  iSw 

CCr^r'CN  /D3LE3/    DUZ  (  5Ci  )  i  DZZ  (bO  i  ) 

CCKNON  /PARAM/    N,NF1,I£X 

CCNNON  /HEAD/    HEC(20) 

CCNMON  /SINGl/    ZM,Ce,CMIN,Cf^AX,  PLTNAX,UM,FG 

CCN'MON  /DMAP/    DKl(lOO) 

CCf^NON  /FLCT6/    RANGEl^) 

CCNMON  /INPUT/    NUNV,NMGD 

CCN.NON  /BGTTO/    DEC  AY,  26  (  20  )  ,U&  {  2C  ) 
C 

€<<<<<<<<«<<<<<<<<«<<<<<<<<«<<<<«<    PROGRAM    BEGINS    >>»>> 
C 
C 


C 
C 

c 


c 
c 
c 


c 
c 
c 


c 

c 
c 


c 
c 
c 


c 
c 
c 


CK=DK1(N) 

***   CHECK  IF  FULL  PLOT  OF  ALL  FCINTS  DESIRED   =*** 

IF  (JUMP-1)  1,2,1 

1  JL^P=NP1/100 

2  JLf'P2  =  JU^P*2 
NLNB=NP1/JUMP 

*^*   PRINT  THE  HEALINGS   ** 

l^PITE  (6,0)  HED,Fe 
WRITE  (6,7)  N,D(<1(^) 
WRITE  (6,8) 

*^*      CONFUTE  THE  RANGE  FACTORS  FCP  PLOT  SIZE   ^^'^ 

PANGEin  =  l.l 

RANGE(2)=-1.1 

R/!NGE(3)=ZN 


**^       WILL  BOTTOM  PROFILE  BE  INCLICEC? 

IF  (ZM-PLTMAX)  4,3,3 

2  NCCCUR=0 

RANGE(4)=G.G 
GC  TO  5 

**=?   DRAW  BOTTOM  FBCFILE   *** 

4  VCCCUR=1 

RANG£(4)=ZN-PLTMAX 

CALL  UTFLCT  (UB , ZE , 20 .RANGE , 1 , MCCCUR ) 
MCDCUR=3 

**♦   DRAW  OCEAN  U(Z)  PROFILE  ^^^ 


** 
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c 
c 


5  CALL  UTPLGT  ( CUZ , CZZ , NUMB , RANGE , JUMP2,M0CCLR) 


FETURN 

FCPMAT 
FCR^'AT 
FCRMAT 
I  • 
ENC 


(•  I', 
CO', 
(  'O'  , 


VS  DEPTH') 


20A4,  '  FREQUENCY  : 
T25,  '^^GDE',  14,  ' 
//,  T23,  'NORMALIZEC 


,  F6.  I,  •  HZ' } 

K  =',  F16.12) 
EIGENFUNCTICN  ' 
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c 
c 


SLBROUTINE    CUTPUT     (CUZ , CK, MODE, NFTS  ) 

Ih'PLICIT    REAL-*8{D) 

CINENSICN  CZZ(NPTS),  DUZCNPTS),  CCCINPTS) 

CCNf^ON  /PHC/  RO,Ra 

CC/^MON  /SINGI/  ZM,C8,CMIN,CPAXtPLTNAX,UM.»FC 

CCFN.GN  /BGTTG/  DEC  AY,  ZB  (  20  )  ,UB  (  20  ) 

CCNMON  /DBLE4/  D S ,D SSCC ,OPVEL , DGVEL 

CCNMON  /HEAD/  HED(20) 

cc^^'0^  /wcrk/  extra(20) 

DATA    NPRINT,NPUNCH/6,7/ 

1  yvRITE    (NPUNCH.IO)     MODE  ,FQ  ,DK,DPVEL  , DGVEL  ,  CS  , DECAY 

DC  2  1=1,20 

2  EXTRA(  I  )=UB(21-n 


c 
c 


c 

c 

c 

c 
c 


WRITE  (NPUNCH,L1}  EXTRA 

IvRITE  (KPUNCHtLl)  CUZ 

WRITE  (NPRINT,13)  MGDE 
RETURN 


ENTRY  PROFIL(CZZ,CCC,NPTS  I 

WRITE  (NP0NCH,7)  HED 

NPTS20=NPTS+20 

WRITE  (NPUNCH,8)  NFTS20, RO, BE, CM  IN t CB, ZM , PLTNAX 

DC  3  1=1, 20 

EXTRACI )=ZN-Ze(21-I) 

WRITE  (NPUNCH,12)  EXTRA 

C2M=DBL£{Zf^} 

DC  4  I=l,IMPTS 
CZZC  I)=D2K'-CZZ{I  ) 

WRITE  (NPUNCH,12)  DZZ 

DC  5  1=1, 2G 
EXTRACI  )=Cfi 


WRITE  (NPUNCH,12)  EXTRA 
WRITE  (NPUNCH,12)  CCC 
WRITE  (NPRINT,SJ 


C 

c 


DC  6  I=1,NFTS 

6  DZZ(  I}=DZM-DZZ(  I) 

RETURN 

7  FCRMAT  (20A4) 

e  FCRNAT  (I10,2F10.5,4F10.2) 

9  FCRMAT  COSCUNO  ViLCCITY  DETAILEC  PROFILE  PUNCHED*) 
END 
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SLtRCUTINE  CZFLQT 
IMPLICIT  REAL=*=8(D) 
C 

C  *  THIS  "SUBRCUTINE  PLCTS  THE  SCLND  VELCCITV 

C  *  PROFILE  CN  THE  PRINTER  PLOT  USING 

C  *  The  N.P.S.  ROUTINE  UTPLCT. 

C  * 

C 

C<<<<<<<<<<«<<<<<<<<<<<<<<<<<<<<<<<<<<<  COMMON  6LCCKS  >>>» 
C 

CCyNiO.M    /DBLI/    DVZ(  501)  ,DH,DH2»DCBSC  ,DR0RB,CK/iF2,DUNAX 
CCNMON    /DeLE2/    D W2 , CCB , DCM IN , DW  , DC  1  ,DC LS G 
CC^MOl\    /DBLE3/    DUZ  (  5C1  )  ,DZZ  (50  I ) 
CCNMGN    /PiiRAM/    N,NP1,IEX 

cc^^;oN  /fead/  hec(20) 

CCf^MON  /UCRK/  WRKL(20) 

CCNMON  /SINGl/  ZM , CB , CM  IN , CMAX, PLTMAX,UM , FG 

CCNf^CN  /PLCT6/  ftANGE(4) 

CCMI^GN  /DSCUND/  DCC(501) 

CCFNQN  /6GTT0/  DEC  AY, ZB ( 20 ) ,UB ( 2C  ) 
C 

C«<<<<<<<<<<<<<<<«<<<<<<<<<<<<<«<<<  PROGRAM  EEGINS  >>>>» 
C 
C      ***   V^RITE  THE  HEADINGS   *** 


C 


C 

c 
c 


t^RITE    (e,5)     nED 
V^RITE     (^»6) 
WFITE     (6,7) 


IF     (PLTMAX.LT.ZM)     PLTMAX=ZM 

Fe=(ZM-FLTMAXi"0.5C-l 
C 

C  -*- DEPTH    INCREMENT    LCCP    BEGINS 

C 

DC    1    1  =  1, 2C 

WFKKI  )=CB 

ZE(n=HE^-^fLOAT(I-l) 
1  CONTINUE 
C 
Q  ^ ,^ DEPTH  LCCF  ENDS 

C 

C      ***   COMPUTE  THE  RANGE  FACTORS   *=** 
C 

RANGE  {  1  )  =  AMAXUCB,CMAX} 

R^NG£(2 )=CMIN 

R/SNGE(3)=2M- 
-C 
C      ***   WILL  BOTTOM  PROFILE  BE  DRAWN?   -** 


IF  (ZM-PLTMAX)  3,2,2 


NCDCUR=0 
RANGE(4)=0.0 
GC  TO  4 


3  MCCCUR=1 

RANGE(4 J=ZN-PLTMAX 
C 

C      =***   DRAW  THE  BOTTOM  C(Z)  PROFILE  *** 
C 

CALL  UTPLOT  ( WRKl  ,  ZE ,20 ,RANGE , 1 , MCDCUR I 

MCCCUR=3 
C 

C      ***   DRAW  THE  OCEAN  CCZ)  PROFILE  *** 
C 
C 

^  CALL  UTPLCT  ( CCC , CZZ ,NP 1 , RANGE , 2 ,MCCCUR ) 
C 

RETURN 
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5  FCPNAT  (MS  20A4 ) 

6  FCKMAT  ( 'O* ,  // ,  T28, 

7  FCRMAT  CO',  //,  Til, 
1  •NETERS  AECVE  BOTTCN 

ENC 


'FLOT  CF  C(Z)   VS   CEFTHM 
•BOTTOM  AT  ZERO.  DEPTH  IN  • 
C(Z)  IN  KETERS/SEC  ) 


I^PLICIT  REAL*8  (A-H,  V-Z ) ,  L0GICAL*4  ($) 
C 

c 
c 
c 

C«<<<<<<<<««<<<<«<<<<<<««<COMMON  BLOCKS  >>>>>>>>>>>>>» 
C 


CCNNCN 
CCNMON 
CCN^'lON 
CC^'MON 
CCNMON 

cc^MON 


/  CI  / 

/  02  / 

/  D3  / 

/  04  / 

/  C5  / 

/  06  / 


CCN^ON  /  07  / 


CCMMGN 
CCf^MOiN 
CCM^ON 
CCi^tVON 

cc^^'QN/ 

CCNI^ON 
CCNiVON 


/  OS  / 

/  010/ 
/  OK  / 
/AWKB/ 
CKi-'AP/ 
/HEAD/ 
/  11  / 


CCMiMON  /  -12  / 


CCr'MON 
CCMMQN 
CCMMGN 
CC^iMO.M 


13 

10 
LI 
PI 


VELG 
ZZZ( 
CC(5 

ORn» 

CM  AX 
H,  CH 
DH2, 
OK,C 
DIES 
ZiMAX 

woce 

F(,V( 

OKI  { 

H£D( 
NUMV 
NP2i 
LOMO 
JHI  , 
JUPP 
NREA 
$GPA 
DPI  , 


(102 
1G21 

Oi,D 
ORB, 
,CMI 
,HPL 
H80T 
Kf^IN 
T,DS 

,css 

2C)  , 

Cl!\T 

100) 

20) 

tNMG 

,NPT 

DE,  I 

JLCW 

ER,  J 

D,NP 

FH,i 

DPIC 


1),  DEPTH  (1021) 
),  DVZ(lOCl) 

EP(50) 

DRORB 

N,CB,CBSQ 

GT,0hhLF,DH2SIX,DH2C3, 

T,DHB 

,DKMAX,DKOLC,CKUf DKL, 

MAX tDE 

tCECAY,DPVEL,DGVEL,G/iy 

NP1,WQCMIN,CVZB 

FCtW,^\2 

,GIFF,OS,DF 


0,NFREC,N^^CCF,N♦NP1» 

S,NCDE,IT  ,NF2 

HNGDE  ,Kl^ELL,IWELL,NWELL. 

BCTT,ISHAF,JTPL, JTPU 
RINT,NPUNCh 
eCTPR,  JCAPCS,$CEL 
2,DPIC4,DPiRT2,D2PI 


C 

C   ' 

C«<<<<<<<<<<<<<«<«<<<<<<<<  PROGRAM  EEGINS  >>>>>>>»>>>>» 

C 

C 

WRITE  (KPRINT,6) 

CALL  INFUT2 

CALL  LINEAR 

IF  ($GRAPH)  CALL  CZPLCT 

IF  ($CEL.OR.$CARDSi  CALL  PRCFIL  ( C EPTH , VELC  ,  NPTS ) 


C 

c- 

C 


DC  5   IFREG=ltNFREw 


FREQUENCY  LCCP  BEGINS 


C 

c- 

c 


FC  =  FQV(  IFREQ) 
W  =  FC>"-D2FI 

wcce=w/CB 
^ccMIN  =  ^^/c^<IN 


CC  1  J=1,NP1 
UCC=W/VELC(J) 
I  DVZ(  J)  =  (^^CCMIN-wCC)*(  l-CC+WOCMINJ 

WCCNPl  =  ^vCC 

CVZB=( ACCyiN-WOCB)*(WCCMIN+WOCB) 

CALL    MAXMI3     (£.5) 


•COMPUTE    VELOCITY    FUNCTION 
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c 
c- 
c 


c 
c 
c- 
c 


c 
c 
c« 

c 


KWELL=1 

LCNCDE=C 


■^'CDE  LOOP  BEGINS 


CC  4  MQDE=1,NM0DF 

C/5LL  SE/iRCh 

C/iLL  ShAPE4 

C^LL  IiNTEGR 

IF  (IGR/iPH)  GO  TO  2 

IF  (MCCE.EC.l)  WRITE  (NPRINT,7)  l-EC 

V^RITE  (KPRINTtS)  N'OCE  ,  CK  ,  DP  VEL  ,  DG  VEL  ,  DSS  ,  CECAY 

GC  TO  3 

2  CALL  MODPLT 

3  IF  ($CEL.OR.$CARDS)  CALL  OUTPUT  (2ZZ,NPTS} 

4  CCNTINUE 

5  CONTINUE 


■FREQUENCY    AND    MODE    LCCFS    END 


STCF 


<<<<<<<<<<<«<<<<<<<<<<<<<<<<<<<<<<<  FORMAT  STATEMENTS  >> 

fc    FCR^AT     CI    NORMQA     :_WKB_METHCD    NOR^'AL    MODE' 
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ELCCK  C 
INPLICI 
CCMVUN 
CCNMQN 
CC^MON 
CCMMOiM 
I 

CCNMON 
CCMMO.N 
D/TA 
CATa 
LATA 
CATA 
D^iTA 
C^TA 
C4TA 
C^TA 
ENC 
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SLEROUTINE  MAXMI3  (*) 

INFLICIT  Rf:/iL*8  (A-H,  V-Z  )  ,  LGGICAL*4  ($) 


LAST  CHANGED  3  AUG  73 


C 
C 

c 

c  

c 

€<<<<<<<<<<<<<<<<<«<<<<<<<<<<<<<<<<<<<  COMMON  BLCCKS  >>>>>> 
C 

CCNMON  /  C5  /  CMAX,CHlN,C6,C.iSQ 

CCFy.ON    /    D?    /  DK,CK(^1K,0KMAX  fDKCLCCKU,  CKLt 

2  DTEST,DSKAXiDE 
CCNNON    /    DK    /  F4.V{20)  ,FC,U,tv2 
CCf-'MON     /AWKB/  DKN  ,  D  I  NT  ,C  I  FF  ,  OS  ,  CN 
CCf^MCN    /HEAD/  HtC(20) 

CCMMON    /    II    /  NUMVTNMOC,NFkEC,NMCCF,N,NPl, 

3  NP21,r.PTS,MGDE,IT,NF2 
CCf^^ON    /     10    /  NREAC,i\PRL\T»NPUNCF 
CCN'MON    /    LL    /  $GRAPr,$BCTPR,iCAHCS,$CEL 
CCMiMON    i    ?l    ^  DPI,DPI02,DPIC4,DPIf<T2,D2PI 

C 

C<<<<<<«<<<<<<<«<<«<<<<<<<<<<<<«<<<  PROGRAM  EEGINS  >>>>> 

C 

c^=o.oco 

DKMN=^/CB 

CKNAX  =  W/CiXIN 

CKN  =  OK-MIN 

CALL    WKfc 

CSMAX=DS/DPI 

NCENAX=  ICIiNT(DSMAX} 

CKCLD=OKMAX 

WRITE  (NPftINT,3)  HED,FQ 

CK^  =  OK.^:AX 

NCDMIN=C 

ksFITE  (NPRira,4) 

CATA  MAXI/'MAXI' 


C 
C- 

c 
c 


c 
c- 


/,     MIM/'MIMV 
MAXI  tOKMAXf  f'GC^'l^ 
WRITE     (NPi^INT,5)     M  IN  I ,  DKM  IN  ,  MODM  AX 


WRITE  (NPRINT,?) 


•CHECK  NUMBER  MODES  PEGU5STED 


NNCCF  =  i\NCC 

iF  (MODNAX-KMQD)  1»2,2 

1  CCNTINUE 


■IF  WE  GET  FERE  WE  MLST  RESET 
THE  NUMBER  CF  MODES  REQUESTED 


KNCOF=MCCMAX 


WRITE 
WRITE 
WRITE 
WRITE 
WRITE 
WRITE 


(NPRINT,6) 
(NFRINT,7) 
(KPRINTtS) 
(NPRINT,?) 
( NPRINTt6) 
(NPRINT, 9) 


NMCDF 


2  IF  (MOCNAX.NE.O)  RETURN 


•IF  WE  GET 
CUT  -  OFF 


FERE  WE  ARE  BELOW 
FREQUENCY 


WRITE  (NPRINT, 6) 

WRITE  (NPRINT,?) 

WRITE  (NPRINT, 8) 

^RITE  (NPRINT, 7) 

WRITE  (NPRINT, 6) 

WRITE  (NPRINT, 10) 

RETURN  1 
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c 

6 
7 

£ 

9 

IC 


FCRMAT 

FCFMAK 

FCRi^AT 

FCFMAT 

FCh^'AT 

FCRMAT 

FCRMAT 

FCBMAT 

1  • 

2  'IjENCV 
END 


•C 

(T33, 

(T33, 

{T33, 

CO', 
I 


HCR12CNTAL  ^.AVE  NUMBER*) 

,14 


F12.8,3>  ,'NCDES 


CIS  lOAe,  4X,  •FRECUENCY  :  ',  Fd.l) 
,//  ,T25,  'BCUNDS  ON 
•t  T22,  A4,  »MUM  K 
i3(  •*'  )  ) 
ih-i,  llXtlH*) 
«-^   CALTIGN   •■*•.«) 
T26,  'MJUES  REQUESTEC  RESET 


(•O* ,/  ,T25, 
CUT-OFF' t//TT20i 
SELECTED'  ) 


'FRECUENCY  RECUESTEO  IS 
JOB  TERMINATED  OR  NEW 


TC  :  ',  I 

EELQW 
FREQ' , 


4) 
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SLERCUTINE  IMEGR 

INPLICIT  REALMS  (A-H,  V-Z),  L0GIC/^L-«*4  ($) 
C      

C      LAST  CHANGED  9  AUG  73 

C      --^^ 

CCNyJN  /  Dl  /  VELD(1021),  DEPTH  (1021) 

CCNMQN  /  C2  /  ZZZ(1C21),  DVZ(lUOl) 

CCPHON  /  D4  /  DROtCRbjDkORB 

CC^MCN    /    C5  /  CMAXtCiy  IN,CE3,CeS& 

CCf^MOK    /    D6  /  H,CH,hFLCT,DHhLF,Dh2SIX,DH2C3, 

1  DHPHPnTTOHP. 

CCN^ON    /    C7    /    DK,CKMINtCKiXAX,CKGLC,CKU,DKL, 

2  OTEST,DSMAX,DE 

CCMMQN  /  DS    /    ZMAX,CSS,CeCAY»DPVEL,DGVEL,GA^ 

CCN/^ON  /  DIO/    k^OCE,l>LCNFl  tWCCr'IN^CVZb 

CC^^'GN  /  Gift    /    FCV  (20)  ,FG,W»l<2 

CCMMON  /  II    /    NUiXV,NMaDTNFR£C,Nf^CCF,N,.MPl  , 

3  NP21, NPTS, MODE, IT ,NP2 


C 
C 


c 
c 

c 
c 


DSS=C.OCO 
C.SSCC  =  O.0D0 


c 
c 

C      <<<<<«<<<<<<<<<<<<<<<<<<<<<<  INTEGRATING  LOOP  BEGINS 

C 

C 

CC  1  1=2, N 
C 

ZZZI2=ZZZ( I) 

IF  (DAES(ZZZI2) .LT.i.OD-30)  GO  TO  1 

ZZZI2=ZZZ12^ZZZI2 

CSS=DSS+ZZZI2 

VELCI  =  VELC(I  } 
2ZZI2=ZZZI2/(V£Lai*VEL0I) 

CSSCC=DSSCC+ZZZi2 

1  CCNTINUE 
C 

C      <<<««<<««<«<<<«<<<<«<«<  INTEGRATING  LCCP  ENDS 

C 

C 

IF  (DABS(ZZZ(NP1 ) )  .LT.l  .OC-30)  GC  TO  3 

A2=DSQRT(CVZ3-DE) 

e  =  ZZZ(NPl)vZZZ(NPiJ*DRaR6--^DRQ/(  2.0DG=^A2) 

2  CSS  =  CRO'=Oh  -CSS+B 
CSSCC=DSSCC^^DRO-OH  +  B/CBSQ 

•  CECAY  =  W*=b/(CB*OK^DSS} 
C 

CPVEL=W/OK 

CGVEL=OSS/(DPVEL^CSSCC) 

RETORN 


3  e=C.ODO 
GC  TO  2 

ENC 
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c 
c 


SLBROUTINE  WK3 

IMPLICIT  RHAL-8  (  A-H  ,  V-2  )  ,  LCCi  C  ^  L'-i-H  {$) 


LAST  CMANG50  :  7  SEPT  73   ASYMP7CTIC  B .  C 


C 

c 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<«<<  CCMMON  ELCCKS  >>»>> 
C 


CCKMQN 
CC/^VUN 


/  D2  / 

/  D4  / 
/  D6  / 


CCMMON  /  D7  / 


CCMMON 
CCM.^  ON 
CC^MON 


/  DIG/ 

/AWKB/ 
/  II  / 


CC^MGN  /  12  / 

CCMMON  /  13  / 
CCMMON  /  LL  / 
CCMMON  /  Pi  / 
CCMMON  /  DTH/ 
CATA  DEXPU  /I 


ZZZ 
ORG 
rii  C 
UH2 
DK, 
DTE 

woe 

OKN 
NUM 
NP2 
LOM 
JHi 
JUP 
$GR 
DPI 
DTH 
.702 


2(1021),  DVZdOOi) 


H,H 
,He 
DKM 
ST, 

,DI 
V,f\ 
liN 
ODE 
tJL 

pr-R 

APh 

,0F 
ETL 


PLC 
OTT 
IN, 
OSM 
CCN 
:JT, 
MOD 
PTS 
,  IH 

ow 
,  ja 

,  $3 
IC2 


kORB 
T,OHhLF,DH2SIX,DH2C3, 

HR 


,DH3 

OK.MAX,OKQLC  ,DKU,OKL 
AX,DE 

,FI,'aOCMIN,CVZB 
DIFF  ,C5,DM 
i,NFf<£C;,NMCCF,N,-MPl, 
., MODE, IT  ,NF2 
IMODEjKwELL  tlVvELL  ,NWELL, 

CTT,ISHAP,JT?L, JTFL 

OTPK,$CARCS,$CEL 
,DPIC4,DFIFT2,D2PI 


DWELL  /6.0/ 


C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<«<<<<<<  PROGRAM  EEGINS  >»» 

C 

CINT=O.CDO 

JTFL=I 

CE=(DKMAX-OKN)*(DKMAX+CKN) 

CKB=DSORT{ DVZ&-Ot) 

$TCP=. FALSE. 

IVvELL=I 

ISIG=-i 

CAPG=O.CCO 

CARGP=C.CDC 


C 

c 
c 

C 


c 

C 
C 
C 


c 
c 

c 


c 

C 


c 
c 
c 


•DEPTH    LGCP    BEGINS 


DC    21    J=1,NP1 
CKZ2=DE-DVZ( J) 


•CHECK  FCR  CONDITIONS 
AND  TYPE  CF  SCLUTIGN 


If  (DKZ2)  1,2,3 

1  IF  (ISIG)  4,5,6 

2  IF  (ISIG)  7,21,9 

3  IF  (ISIG)  10,19,20 


IMAGINARY  REGION 


CAPGP=DSQRT(-DKZ2) 
C4PG=CARG+CARGP 
GC  TO  21 

ISIG=-1 

GC  TO  4 

C/!RGP  =  DSCRT(-DKZ2) 

DX=CKZ/ (CKZ+DARGP) 

CINT=DINT.+  CKZ='(DX-i.OD0)-0.5DO^CF 

ISIG=-1 

JTFL=J-1 

CARG=(  2.0DC-DX)'^DARGP*C.5D0 

GC  TO  21 


ZERO  VERTICAL  WAVE  NUMBER 
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c 

c 

c 


c 
c 
c 
c 

c 
c 
c 


10 


u 


12 


12 

1^ 


it 


1? 

18 


19 


2C 


/Lc 


23 


2A 


ie/5CK=-l 
ISIG=0 
GC    TO    ^1 

IE/^CK=I 
GC    TO    8 


•REAL    REGION 


IS 
IF 
JT 
CK 
CX 

IF 
CT 
CI 
ST 
GC 
CT 
IF 
IF 

IV. 
JT 
CT 
CI 
GC 

e  = 

CE 
Cc 
f^C 
EC 
DI 
IF 
Ci 
GC 
CI 
CI 
C/i 
GC 


IG=l 

($TC 
PC  =  J 
Z  =  DSU 
=  CARG 
KG=DA 

($Tu 
HETU  = 
NT  =  DT 
CP=.T 

TO  1 
1  =  DM0 

(DAR 

(IWE 
PG=C. 
ELL=I 
Pb=J 
hETU  = 
NT  =  DT 

TO  L 
C3I.\( 

ccas( 

L  =  OEX 
fa  =  DE 
=  CEL^'= 
=  DEL- 
NTP=D 

(DIN 
NTP  =  D 

TO  L 
NT  =  DI 
KT  =  OI 
RG=0. 

TO  2 


F)  GC  TO  11 


RT( 
F/( 

kG-» 

P) 

CAT 

FET 

PUE 

e 

C(D 
G=^D 
LL. 

ceo 

UEL 


DKZ2i 

LKZ+OARGF) 

DARGP^^(DX-I.CC0)*C.5DC 
GO  TO  13 

AN{DTANF(CARG*CH)  ) 
U 


INT,D2PI  ) 
h-CWELL  ) 
EG.KWELL)  GO 

L  +  1 


TO  ZQ 


CP 

FE 

£ 

DT 

CT 

P( 

XP 

(A 

(A 

in 

TP 
IN 
fc 

^T 

NT 
CC 
1 


104 
TU 


11 
1  ) 

DA 
(- 

+  e 

+  6 
T- 
.G 
TP 


RG^DH) 

DARG --Dh) 

)  +  Ccf^L+(A-B) 

DT1+DATAN2( AC,BO) 
E.DINT-CTtST)  GO  TO 
+DPI02 


17 


+  {2.CO0-CX)-DKZ'^C.5D0*=DH 
0 


CKZ=DSQRT(CKZ2) 

ISIG=+i 

CINT=DINT+DKZ^DH 

IF  ($TCF)  GC  TO  13 

JTPC=J 

GC  TO  12 

CKZ=DSCRT(DKZ2) 
CINT=DINT+OKZ'»DH 


21  CONTINUE 


IF  (DKZ2)  22,27,25 


CEPTH  LCCF  ENCS 


•DECAYING  SCLLTICN 


CKS=DARGP 

eARG=(DA^G-C.5D0'i«CARGP  )=feDH*2.0D0 

IF  (CAKC-OEXPUi  23,24,24 

C^=(CKS-i')R0R6"-0KB)*CEXP(-DARG}/  (  CKS  +  DRORE*C  KB  ) 

CTFETL=CATAN((1.0CC+D2)/tl.GD0-D2)) 

GC  TO  27 

CThETL=CPI04 

GC  TO  27 
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c 
c 


JTFL=NP1 

IF     (DKB.tC.C.ODO)     GO    TO    26 
CThETL  =  CATAN(DKZ/(CKB'f=CRGRB)) 
GC    TO    21 


BOTTOM    REFLECTED 


2e    CTI-ETL  =  CFIC2 


■COMPUTE    CIFFERENCE 


21 


C£  =  CINT+QTI-ETL 
CINT=DINT-DTHETU 
CIFF  =  DS-DM-*DPI 
N^;ELL  =  IWELL 
FETURN 


6     ItNELL=IV^ELL+l 
•GC    TO    24 


ENC 
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SLEROUTINE  SEARCH 

I^FLICIT  i^E^L    ^8  (A-H,  V-Z),  L0GICAL*4  1$) 


LAST  CHANGE  :  15  AUGUST  1973 


C 

c 
c 

c 

C<<<<<<<<«<<<<<<<<«<<<<<<<<<<<<<<<<<<  COMMON  BLCCKS  >>»» 
C 


CCMMON  /  C7  / 


CCNi'^ON  / 
CC^NGN/CK 
CC^MCN  /A 
CCMMON  / 


D6    / 

iVAP/ 
WKB/ 
I  1    / 


CCNI/QN    /    12    / 


CCNMQN 
CC^MON 
CC/^HON 

cc^^<CN 

CCMiMGN 


13 

IG 
LI 

PI 


DTH/ 


DK,C 
DTES 
DIFF 
DKK 
DKN, 
IMUMV 
NP21 
LOMC 
JHI  , 
JUPP 
NREA 
$GRA 
DPIt 
DTHE 


K^'IN 

T,DS 

I 

100) 

C  I N  T 

,iNHG 

,NPT 

CE,I 

J  LOW 

ER,J 

DtNP 

PH,$ 

CPIO 

TU 


tDKMAX,DKOLC,CKUiOKL, 
MAXtuE 


t  C I FF , CS  »  C^ 
D,NFREG,NNCCF,N,NP1 , 
S,MUDEt IT,NP2 
HNGDE  ,KWELL,IWELL,KWELL, 

eCTT,  ISHAP,JTPL,JTPIJ 
RINT,i\PUfMCF 
BCTPR,$CARCSi$CEL 
2tUPIC4,0FIRT2,02PI 


C 
C<<< 

c 
c 
c 
c 


CATA  ITNAX  /30/ 
<<<<<<<«<<<<<<«<<<<<<<<<«<«<<<<  PRGGRAV  BEGINS  >>>» 


C 
C 
C 


C 
C 


c 
c 
c 
c 


FIX=1.400 

IF  (imol;e.gt.u  go  TG  1 

KWELL=1 

$LCh=. FALSE. 

$hl=. FALSE. 

CKCLD=DKMAX 

JLPPER=  1 

Jhl  =  l 

JLCW=NP1 

IFMCDE=C 

LCNCD£=C 

CKCEN=CKGLD 

CKC£N2=CKCEN 

jeCTT=NPl 


INITIAL  SETTINGS 


■SET    UP    FOR 


DM=DFLGAT(MODE-LaMODE-IHMGDE) 

IT  =  0 

IF  (SLOVn}  GC  TO  11 

CKGUES=(CKCLD-OK.M  IN) /(DSMAX-DM+ 1.000) 

CKL=CKCLC-CKGUES*FIX 

IF  (OKL.LT.DKMIN)  CKL=CKMIN 

CJ<U  =  CKQLD 

FL=-OPI 

CKN=DKL 

CALL  WKB 

IF  (JTPL.LT.JHI)  GC  TC  5 

IF  (OIFF.LT.O.ODO)  GO  TC  4 

FL=CIFF 

JHI=JTPL 

JLCW=JTPL 

GC  TO  16 


GUESS 


4  FIX=FIX+C.4D0 
GC  TO  2 


'UPPER  CHANNEL  SET-UP 
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c 
c 

c 


c 
c 
c 

c 
c 
c 


c 
c 
c 


5  IF     (DS-CFLCATdHMCCE  +  l  }*DPI  ) 

6  KU£LL=2 
GC    TC    3 

32     $Fl=.Tf<Uf:". 
JLFPbR=JTPL 
IF^CDt=Ih^'CDE-»■l 
C/'  =  DFIO/iT(  iHMQOc  J 

7  CKN=(CKU  +  DKL)'^0.5D0 
6    Cy^LL    WKB 

IF  (JTPL.LT.JUPPER)  GO  TO  9 
CKU=DKN 
GC  TO  7 
S     IF  (DIFF.LT.C.ODC)  GO  TO  10 
CKN=(CK0+DKN)=^0.5C0 
GC  TO  8 
IC  CKU=DKN 
FL=DIFF 
GC  TO  16 


6,32,32 


11  CKL=DK 

CN  =  CFLC/iT(LCMODE) 

CI<L  =  DKCEN2 

CKN=(DM  +  DKL)*0.5CC 

12  C/iLL    WKE 

IF    ( JTPL.GT.JfiOTT)     GO    TC    13 

DKN=(  CKt-»-DKNJ^'0.5DC 

GC    TO    12 

FL=DIFF 

CKlj  =  DKN 

IF     (FL.GT.O.OCO)     GC    TO    16 

CKL=DKL-.21iDC-(L)KyAX-DKMi\)/DSM/iX 

IF     (DKL  .LT^DKMIN)     CKL=OKHIN 

CI<N=OKL 

C/iLL    l^KB 

IF     (JTPO.LT.JLOW)     GO    TO     15 

FL=CIFF 

GC    TO    1^ 

K^sELL=i 

$LCW=.F^LSE. 

LCNCDE=LGMC0E-1 

GC  TO  1 


13 

14 


15 


17 


le 

IS 


LOWER  CHANNEL  SET-LP 


COMPUTE  NEXT  DKN 


le    CKN  =  DKL  +  FL''(CKU-CKL)/(FL-FU) 


'CALL  WKB 


CALL  rJKB 

IT=IT+I 

IF  ($LOV..OR.$HI  )  GC  TC  18 

IF  (JTPt'.LT.JUPPER  )  CIFF  =  DIFF-DFLGAT(IHMCCE)*DPI 

IF  (JTPL.GT.JBCTT)  DI FF=D I FF-DFLOAT ( LOMOC E J *CP I 


IF  (OIFF)  19,20,21 

CKl;  =  DKN 
FL=DIFF 
GC  TO  Z2 


•CHECK  WHETHER  ABOVE  CR  BELOk* 


2C    WRITE     (NPRINT,31)     ^'CDE 
GC    TO    25 

21    CKL=DKN 
FL=DIFF 


•TEST    THE    RESULT 
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c 
c 
c 


22 


23 


24 


25 


2t 


C 

c 
c 


21 


28 


29 


C 
C 
C 


3C 


CINC  =  D/i3S(CKGLD-DKN)/DKN 

CKCLD=DKN 

IF  iDAbS(CIFF)-DTEST)  25t25,23 


IF  (MOC{IT,IO) .EC.9)  GC  TO 
IF  ( IT.LE.ITMAX)  GC  TC  16 
IF  (OINC.LT. I .00-15)  GC  TO 
CKN=( CKL+DKL)«J.5C0 
GC  TO  17 


d^ 


•THE    NCOe    H/^S    BEEN    FCUNC 


CKUHOCfc  )  =  CKN 

CK=CKN 

CIFFI=DIFF 

IF    (5hl  )    GC    TO    26 

IF    ($LOW>     GO    TO    2  7 


CKCEN2=CKCEN 

CKCEN=CKK 

JFJ=JTFL 

JLCW=JTPL 

IF     (JTPU.LT.JUPPER) 

IF     (KWELL.LT.NWELL) 

KtsELL=I 

IF     (JTPL.GT.JEOTT)     LONCCE=0 

ISFAP=I 

RETURN 


•CENTER    WELL 


IHf^CCE  =  0 
GC    TO    28 


•LPPER    WELL 


CI'CLC  =  CKCEN 
JLPPER=JHI 
$hl=. FALSE. 
K^^ELL  =  2 
ISFAP=2 
RETURN 


•LOWER    WELL 


CI<CLD  =  CKCEN 
KWELL=l 
$LCW=. FALSE. 
ISFAP=4 
RETURN 


K^ELL=KWELL+1 

JShAP=3 

J1=JTPU 

CS1=DS 

J2=JTPL 

CALL  WKE 

LCCIFF  =  CS-CFLCAT(LC^^CCE  +  l  )'pCPI 

IF  (LOCIFF+IO.ODC^DTEST)  29,30,3C 

KWELL=1 

JTPL=J2 

CS=CSl 

JTFU=J1 

RETURN 


•CHECK  FOR  NODE  IN  LCWER  DUCT 


■SET  UP  FOR  LOWER  CHANNEL 


LC^CDE=LCMCDE+l 

JECTT=JTFU 

FL=LGDlFF 

ILCW=.TFUE. 

JTPU=J1 

CS=DS1 

JTFL=J2 
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RETURN 

31  FCFN/^T  (•OI^CDE  •  ,  14 ,  •  FCUND  BY  CIPECT  HIT.') 
ENC 
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SLEROUTINE  ShAPEA 

IMPLICIT  REAL*i=8  (A-H,V-Z),  LCGICAL*^  ($) 
C 

C«<<<<<<«<<<<<<<<<<<<<<<<<<<<«<<<<<<    COMMON    ELCCKS    >>>>» 
C 


CCf^i^GN 
CCKMCN 
CCNNON 
CCNKON 


/  D2    / 

/  C4    / 

/  C5    / 

/  ce   / 


CCi^MGN    /    D7    / 


CC^.MON 
CCNMQN 
CCf^MON 
CCNMGN 
CCNMQN 


/  C9  / 
/  010/ 
/    DW    / 

/AkvKB/ 
/    11    / 


CCNNCN    /    12    / 


CCNFGN 
CCNMON 

CCNKON 


/  13  / 

/  IG  / 

/  LI  / 

/  OTH/ 


ZZZ 

ORG 
CMA 
H,  C 
DH2 
OK, 
DTE 
ZMA 
UGC 
FQV 
DKN 
NUf 
NP2 
LGM 
JHI 
JUP 
NRt 
$GR 
OTH 


(  1C2U 

tCRE.D 
XtCf^lN 
H,hPLC 
,HBGTT 
CKMN  , 
ST,DSM 
X  ,CSS, 
B  ,hCCi\ 
(20, F 
,CINT, 
V,N^GD 
1,NPTS 
CCE,  Ih 
,  JLQW 
PER, JB 
AD,NPR 
APh, $6 
ETU 


,  DVZ(LOCl) 

RORb 

,CB,CESQ 

T,DHHLF,DH2SIX,0H2Ca  , 

,0H3 

DKNAX,DKCLi:,CKU,CKL, 

AX,DE 

C£CAY,CPVEL,DGVEL,G/iN  . 

PI  ,WQCN'JK,CVZB 

Q  ,  W  ,  W  2 

CIFF,CS,  Ct* 

,NFR£C,NNCCF,N,NP1, 

,MUCE, IT ,KP2 

NODE,  KW  ELL,  I  WELL,  NVv  ELL, 

GTT,  ISHAP,JTPL, JTPL 

IiNT,i\FUNCF 
GTPR,$CARDS,$CEL 


C 

C«<<<<<<<<<<<<<<<<<<<<<<<<<<«<<<<<<<<  DATA  >>>>>»>>>>»» 

C 

C;iTA  CUPPER  /9.99C  60/ 

<<<<<<<<<<<<«««<<<<<<<<<<«<<<<<<  PROGRAN  BEGINS  >>>>> 


C 

c<< 

c 


CSTAR 
CE  =  (D 
IF  (I 
IF  (I 
$GC  =  . 
ISTAR 
ISTCP 
GC  TO 
ITCP= 
GC  TO 
ITCP= 
GC  TO 
ITCF  = 
GC  TO 
ISTCP 
GC  TO 
ISTCP 
CCNTI 
AVKZ  = 
ZJST  = 
F/iCTG 
ZPJST 
ZJ  =  ZJ 
2FJ  =  Z 


T  = 

K,v 
SF 
SF 
FA 
T  = 
=  N 

( 
1 

4 
JL 

4 
JU 

( 
=  J 

7 
=  J 
NU 
DI 
DS 
R  = 
=  A 
ST 
PJ 


l.CCO 

/iX-CK) -(DKN-^X  +  CK) 

AP..EQ.3.A^C.  IhiVCCe.GE.l  )  ISFAP  =  6 

/^P..  EQ.l.ANC.IFMGDE.GE.l)  ISFAP=5 

LSE. 

JTPL 

PI 

1,1,1,2,3,3],  ISFAP 


Cl^  +  1 

FFER 

7,5  ,6,7,7,6)  ,  IShAP 

FI-1 

bCTT-1 

E 

NT/(CFLQAT{ JTPL-JTPU )*Ch) 

Ii\(DThETU)^'^CSTART 

ZJSJ^^iUll  ISTART)-DVZ(  1ST  ART- 1 )  )  /  (  A  VKZ*4.  DO=»'DH) 

VKZ*DCCS  (DTh  ETU  )-^OST  ART -FACTOR 

ST 


C 

c 
c 
c 


FJ=DE-DVZ( ISTARTJ 

JSTART=ISTART 
JSTCF=ISTCP 

CEK0M=i.0DC-DH2SIX*{DE-DVZ( JSTART-l)  ) 
ZJ;a=(  (  1.0D0-DH2D3^FJ)>^ZJ-DH^-ZFJ)/CEN0M 


ZNAX=G.CD.O 

CC  10  J=JSTART, JSTCP 
ZZZ( J)=ZJ 
FJ=DE-DVZ{ J) 

ZJF1=(2.0D0-DF2*FJ)«ZJ-ZJMI 
ZjM  =  ZJ 


•DEPTH    LOOP    BEGINS 
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c 
c 
c 
c 


c 
c 
c 
c 


c 

c 
c 


/S2J  =  CABS(ZJ) 

IF     (J-JTPL)    9,9,8 

e    IF     (AZJ.LT.ZFiAX*CTEST)    GO    TC    11 
IF    (DAbS(ZJPl)  .GT.AZ-J)     GQ    TC     11 

S    IF     (AZJ.GT.DUPPER)     Gu    TC    29 
IF     (AZJ  .GT.ZMAX)     Z^AX^AZJ 
IC    ZJ=ZJP1 


IF  (ISTCP.EC.NPl)  GO  TO  16 

J5TART=JSTCP+1 

GC  TO  U 

11  JSTART=(J+jTPL)/2 

12  ZJ=ZZZ( JSTART) 

CC  15  J=JSTART, JSTCP 
2ZZ( v)=ZJ 

CKZ  =  DSWf<T(DA6SlDVZtJ)-DE)) 

2J  =  ZJ--^DEXP('-CKZ*CF) 

IF  (DAoS(ZJ)  .LT.Zf^'AX^CTEST) 


•CEPTH  LOOP  ENDS 


12  CONTINUE 


G3  TC  14 


lA 


15 


16 


17 

18 

19 

2C 


21 


^^ 


IF  (I5TCF.EC 
J=JSTOP 
JSTART=J+1 
$GC=.TRLE. 


NPl)  GC  TC  16 


CC  15  J=J.STAFT,NPTS 
ZZZt J)=C.ODO 


CCMFUTE  UPPER  DECAYING  SOLUTION 


IF     (ISTART.EQ.ITQP)    GO    TO    24 

2J=ZJST 

ZPJ=-ZPJST 

FJ  =  DE-uVZ(  ISTART) 

ZjM  =  ZJ-Lh-ZPJ-DH2=«'FJ*0.5D0 

ISTAR1=ISTART+1 

JENC=IbTARl-ITOP 

2^XDT  =  ZNAX=?«DTEST 

CC  18  J=1,JENC 

INCEX=ISTAR1'-J 

ZZZ(IND£X)=ZJ 

FJ  =  DE-CVZ{  INDEX) 

2JF1  =  (2.0D0-DH2=^FJ)=^2J-ZJM1 

ZJM  =  ZJ 

/ZJ=CAbS( ZJ) 

IF     (AZJ-Z.f^XDT)    20,20,17 

IF     (DADS(i;jPi)-AZJ)     18,18,22 

ZJ=ZJP1 

IF  (ITOP.EG.l)  GQ  TO  24 
KSTGF=INDEX-1 

CC  21  K=1,KSTCP 
2ZZ{K)=C.0D0 

GC  TO  24 
JSTART=(J-H)/2 

2c=ZZZ( ISTaRI-JSTAPT J 

DC    23    J=JSTART, JENC 

INC£X=iSTARl-J 

IIH  INDEX)=ZJ 

CKZ  =  DS(.RT  (CABS(DVZ(  INCEX)-DE)  } 
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c 
c 
c 
c 


c 
c 


22 


ZJ  =  ZJ=^CEXP{-DKZ=«DH) 

IF  (CABS(ZJ)  .LT  .Zf^XCT)  GO  TC  20 

CCMINUE 


GC  TC  19 


—NORMALIZE  THE  FUNCTICN 


2^  CCCEF=L  .ODO/ZMAX 


25 


CC  25  K=ltNPI 
ZZZ(K)=ZZZ{K)^DCaEF 


' — COMPUTE  THE  BOTTCN  SCLUTICN 


26 


27 

28 


29 

2C 

31 


IF  (  .NCT.iBOTPR)  f-ETURN 
EXFCN=DEXP('CFB»DSGRT(CVZB-CE) ) 
ZZZeaT=ZZZiNPl}^ORCBt 

CC  26  J=NP2»NP21 

ZZZ( J}=ZZZdCT 

IF  (DAbS{ZZZtiGT)  .LT.l.CD-60)  GO  TC  27 

Z2ZEuT=ZZZECT':=EXPGN 

RETURN 


CC  28  I=J,NP2l 
ZZZ(n  =  C.Q 

RETURN 


IF  (DSTART.LT.l.CD-60)  GO  TC  30 

CSTART=CSTART=<=l.OC-lO 

GC  TO  7 

WRITE  (NPRINT,31) 

J=l 

GC  TO  27 


FCRMAT  CI  EXPONENT  PPCTECTION  EXCEDED  IN  SHAPE4',//, 
1      •   RETURNED  TC  CALLING  PROGRAM  INCQNFLETEM 
ENC 
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SLERCUTINE  OUTPUT  (CUZ,NPTS) 

INFLICIT    REiiL*8     (A-h,     V-Z )  ,     LCJGIC/M^4    ($) 

C 

CINBNSICN  DZZINPTS),  DUZ(NPTS),  CCC(NPTS) 

C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  CCMMCN  ELCCKS  >>>>>> 
CCN/^GN  /  04  /  DRQ,CRB,DR0R3 
CCN^ON  /  C5  /  CMAX,CMI,\i,CBtCBSG 
CC^NGN  /  06  /  H,CH,HPLCT  ,DHHLF,DH2SIXiDH2C3, 

1  DH2,HP0TT,DHB 

cc^^'.c^l  /  o?  /   CK,CKMK,L;K;^Ax,OKCLC,CKU,DKLf 

2  DTEST,OS/-*AX,0E 

CCKMCN    /    09    /    ZMAX,CSS,LtCAY,DPVEL,DGVEL,G/i^ 
CCNf^ON    /    DU    /    FQVdC)  ,FC,WtW2 
CCNMGN    /AWK6/    DKN , 0  INT , C IFF  ,DS , DN 

cc^^<ON  /fead/  hec(2G) 

CCFNCN  /  II    /  NUiMV,Kf^CC,NFREC,NyCCF,N,NPl, 

3  NP21,NPTS»iMGDE,  IT,NF2 
CCNMGN  /  IG  /  NR£ACtNPRINT,NPUNCF 
CCNMON  /  Li  /  $GRAPH,$eCTPR,$GAPCS,$CEL 

C 

£<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  PROGRAM  BEGINS  >»» 

C 

C      /^CDE  OUTPLT  SEGTICN 

C 

1  WRITE  (^PU^CH,6)  NGCE,FG,CKN,CPV£L,DGVEL, CSS, DECAY 
WRITE  (NPUNCH,7)  DLZ 

C 

IF  (  .NGT.SGRAPH)  GC  TO  2 

WRITE  (NPRlNTjS)  f'CCE 

RETURN 
C 

2  WRITE  (NPRINTilOJ 
RETURN 

C 

Q  SCUND  PRCFILE  OUTPUT 

C 

E^TPY  PFCFIL  (DZZ,  CCC,  NPTS) 

WRITE  (NPUNGH,3)  HEC 

WRITE  (NFUNCh,^)  NPTS , ORG , CRB, GN  I^ , CB, H, HPLCT 

WRITE  (NFLNCH,6>  CZZ 

WRITE  (NFUNCH.B)  CCC 

WRITE  (^FRI^T,5) 
RETURN 
C 

c 

C«<<<<<<««<<<<«<«<<<<<<<<<<<<«<<<  FORMAT  STATEMENTS  » 
C 

3  FCR^AT  (10A8) 

^  FCRMAT  niCt2Fl0.5,4FL0.2) 

5  FCPI^AT  (•CSCUND  VELOCITY  DETAILEt  PROFILE  PINCHED') 

6  FCRMaT  (I5,F5.0,G2G.12,2F10.2,2G15.8) 

7  FCRNAT  (1CF8.5) 

8  FCFNAT  (10F8.1) 

S  FCRMAT  CO   GLTPUT  FOR  f'ODE*,  14, •  HAS  BEEN  PUNCHEC) 

IC  FCRMAT  (•+•,  T118,  'OLTPLT  PUNCHEC*) 
ENC 
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SLERGUTINE  INPUT2 

1^FLICIT  REALMS  (^-H,V-Z},  LOGICAL-4  ($) 
C 

C        :- 

C  L/5ST    Ch/iNGElD    3    AUGbST     1973 

Q  

c 

C<<<<<<<<<<<<<<<<<<«<<<<<<<<<<<<<<<<<<  CCMMON  BLCCKS  >>»>> 
C 

CCNFGN  /  C3  /  CC(5C)  tCEP(50} 

CC^^'C^  /  D4  /  DRC  tCSBfChOKb 

CCNMON'  /  C5  /  Ci^'.A>(,C^IN,C6,GBSQ 

CCKNON  /  Dfc  /  H,Ch,HFLCT,DhhLF,Ch2SlXTDh2C3, 

1  DH2,hEC7TtOHB 

CC^MG'N  /  D7  /  DK,CKyiN,CKMAX,OKCLC,CKUfOKLt 

2  DTtST,OSMAX,De 
CC^^f^.ON    /    DU    /    FQV(2C),FC,l-<,^2 
CCM^CN    /HEAD/    HEC(20) 

CCMMQN    /    IL    /    NUMV,NM0D,NFREG,NMCDF,N,NP1, 

3  NP21,NPTStMGDE,  iT,l\F2 
CCNMCN  /  IC  /  NRtAC,NPRINT,NPUMCF 
CCNMCN  /  LI  /  $GRAPH,$BGTPk»3.GARCS,$CEL 

C 

C<<<<<<<<«<<<<<«<«<<<<<<<<««<«<<<  CATA  >>>>>>>>>>>>>» 

C 

CATA  CC^VRT  / 1 . OOOCCOOCCOCOOGOO CO/ 

c 

C«<<<<<<««<<<<<<«<<<<<<<<<<<<<«<<<    PROGRAM    BEGINS    >>>» 

C 

C 

C      -— ' ^ READ  IN  THE  HEADING 

C  AND  RUN  FARAMcTERS 

C 

RE^D  (NREAC,9)  HED 

ksRITE  (NPRINT,13J  FED 

READ  (NREAD,10}  NUNV,  NiVCD  ,  I  EX  ,N  ,  ICLT  ,  IDEV 
C 
C      . . '^ TEST  THE  RUN  PARAMETERS 

C 

IF  {N.GT.IOOG)  N=1G00 

IF  (NUMV.GT.5C)  ^L^'\/  =  50 

IF  (NMOC.GT.LOO}  NNCD=  100 

IEX=IABS(IEX) 

IF  (IEX.GT.14)  IEX=14 

DTEST=LC.ODO-^v(-lEX} 

IF  (IQGT.LE.2)  SGPA FH=  .TRUE  . 

IF  ({ IOLT.Ew.2)  .OR.  (ICUT.EQ.3))  $CEL=.TRLE. 

IF  ((  lOUT.EQ. i)  .OR.  (  IGGT.EQ.4)  )  $CAkDS= . TR L E  . 

IF  (lOEV.NE.O)  NPLNGh=iCEV 

NF1=N+1 

NF21=N+21 

BEAD  (NREADtll)  H , CB , DPC  ,DRB ,HPLCT 
C 

C  ■■ READ    IN    THE    FREQUENCIES 

C 

RE/5C  (NREADtLO)  NFPEw 

IF  (NFREG  .GT.20)  NFREG  =  20 

READ  (NRcAOtLl)  ( Fu V ( I )  ,1 =1  ,NFR EG) 

FG  =  FQV(  1) 
C 

C      . CCNVERT  AND  DISPLAY 

C  THE  RUN  FARAyETERS 


C 


IF  (NUMV)  1,2,2 

CCNVRT=C.3C48 

NLNV=-i\LNV 

h  =  F-'CCNVRT 

CE=Cd^CCnVST 

FPLCT  =  hPLCT='CONVRT 

WRITE     (NFRINT,14) 

l^RITE     (KPRINT,15)     (  F3V  (  I)  ,  I  =  1 ,  NFR  EG) 

l^FITE     {NPRINT,16)     NNGD,IEX 
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c 
c 
c 
c 
c 


URITE  (NPRINT,17) 

WFITE  (^PHI^T,la} 

^^RITL-  (^PRIl\TtL9) 

IF  (IGRAPh)  WRITE 

IF  (SCARCS)  ^;PITE 

IF  (iCEL)  WRITE  (NPRINT,20)  IDEV 

WRITE  (NPRINT,23)  HED 


H,HPLCT 

ce  ,N 

(NPRINT,22) 
(i\FPlNTt21J 


•READ  IN  THE  DEPTHS 
AhD    SOUND  VELOCITIES 


C 

c 
c 

c 


CC  2 

RE/C     (NRtAClU     CEFI,CCI 

D£PI  =  DEFr-CCNVRT 

CC  I  =  CCI^CONVRT 

CC  (I)=CCI 

CEP{  n^DEPI 

WRITE  {NPKINT,12)  DEPI,CCI 

CCNTINUE 


CF  =  h/DFLCAT(f>l  ) 

DPCR6=LhC/CRE 

CF2  =  DH-'CH 

CeSC  =  CB''CB 

Dh2S  IX  =  Ch 2 -^0.1666  0666 6 6£66 7 

CF2D2=DH2SIX^2.CDC 

CFFLF  =  t(-^^  0.500 

IF  (DEPI-H)  4,6»5 

NLNV=f\U^V+l 

CC(NUM\/)-{h-0EPI  )  ^'0.016  500  +  001 

l;ep(nomv)=h 


COMPUTE  SCNE  CONSTANTS 


6  NFTS=NP2L 
IF  IHPLCT-h)  7,24,8 

7  hPLOT=H 

24  5eCTPR=  .FALSE. 
NFTS=NP1 
RETURN 

8  HECTT=hPLGT-H 
CFe=HBGTT^'0.5D-l 
RETURN 

C 

C 

C«<<<<<<<<<<<<<«<<<<<<<<<<<<<<<<«<<<  FORMAT  STATEMENTS  >> 

C 


9 
IC 
LI 
12 
13 
14 
15 
16 


17  F 
I 

18  F 
1 

19  F 


20 
21 
22 
23 


1 


CPNAT 
CRMAT 
CRMAT 
CR^AT 
CRMAT 
CRMAT 
CRMAT 
CRMAT 

T24 
CRMAT 

T22 
CRMAT 
•  METE 
CRMAT 
'NUiMBE 
CRMAT 
CRMAT 
CRMAT 
CRMAT 

•DE 
NC 


(10A3 
(5II0 
(dFlO 


0-, 
1', 

0«, 
0', 
u', 
ITE 

C  , 
.«MA 
0«- 


RS/SEC 

( 

R 

( 

( 

( 

( 

PTH 


0', 
,// 
OS 
0', 


) 

,12) 
.3) 

T25 
lOA 
T34 
T2  3 
T24 

RATI 
T27 

X  DE 
T18 
•  ,/ 
T22 

,T42 
T42 
T42 
T42 
LOa 

METE 


,  Fc.l,  T57,  F6.1J 

8) 

,  'RU 
,  'FR 

,  "MO 


ON  LI 

,  'dC 

PTh  p 

•BC 
,Ti7 

•CU 

•GRC 
•GU 
'  FU 
•SO 

// 


FS* 


N  PAR 
ECU  EN 
DES  R 
MIT  : 
TTCM 
LGTT5 
TTCM 

,  •  ^  u  M 

TPLTS 

UP  AN 
TPUT 
NCFEC 
Lf^D  V 
,  T33 
T44, 


AMETE 
CIES 

EgUES 

OEPTH 
D  :•  , 
SOUND 
BER  C 

ftcOU 
D  PHA 
ON  FI 

CARC 
ELGCI 
,  'Vc 
•SCUN 


1, 
,  //, 
MM 

•  VETE 


:  •  .(T42,F6.1,  •  HZ* ,/)  ) 
TEC  :•  ,  14 
(  -M  12, 
:•,  F3.  1, 

f8.l,  •  NETEFS*  ) 

VELGCITY  :«,F8 
F  GRID  POINTS  ;' 
ESTED  :  HCRIZ  WAVE 
SE  VELCCIT  lES'  ,// J 
LE  FT' ,12,  'FCCl') 

CUTPUT'  ) 
TV  AND  MCCE  GRAPHS  '  ) 
LCCITY  PkCFILE' ,//  ,T22, 
C  VELOCITY,  K/SEC,  /) 


RS',//, 

1, 
,15) 
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SLBROUTIKE  MODPLT 

INPLICIT  REAL-^a  (A-H,  V-Z )  ,  LOGICAL-^4  ($) 

CINENSICN  RANGE  (4) 
C 

c  => 

C      *    THIS  SUSROUTIKE  PLCTS  THE  MCCE  SHAPE  CN  THE  PRINTE 
C      *    UTILIZING  THE  STANCARD  (AT  NPS)  ROUTIiNE  UTPLCT 
C      '^ 

c  * 


c 


c 
c 


CCNNON  /  Ci  /  VeLC(1021J,  DEPTH  (1021) 

CCf^MON  /  D2  /  ZZZ(iC21),    DVZ(L0CIJ 

CCM^uN  /  C5  /  CMAX,CMlN,CB,CbSg 

CC^'^'0^  /  D6  /  H,Ch  ,HFLCT  ,DHHLF,CH2SIX»DH2C3, 

1  0H2?HB0TT,DHB 

CCNMON  /  D7  /  DK,i:Ky,IN,CKi^AX,DKCLC,CKUf  DKL, 

2  DTEST,DS^^AX,DE 

cc^MQN  /  ca  /  diffi 

CC^"^'ON  /  D9  /  ZMAX, CSS, DECAY, DPVELjDGVEL, CAN 

CCI^NON/DKMAP/  DKUICO) 

CCNMGiN    /    OW    /  FQV(20)  ,FQ,W,W2 

CC/^MON  /AWKB/  DKN' ,  D  INT  ,  CI  FF  ,  CS  ,  CN 

CCFNON  /HEAD/  HE0(2C) 

CCNMQN  /  II  /  NUMV,NNCC,NFREC;,NNCCF,N,NP1, 

3  NP21,NPTS jMGDE,IT.NF2 
CCMMON  /  10  /  NREAD,NPRfNT,NPUNCH 

WRITE  (NPRINT,!)  HEC,FC 

WRITE  (NPRINT,2)  f^.GCE,DK 

WRITE  (NFRINT,3) 

CALL  UTPLGT  ( ZZZ , DEPTH , NPTS , RANGE  ,2 , 0 , PBCTT ) 

WRITE  (i\PRINT,4)  CPVEL  ,  CSS  ,  DEC  AY  ,  CGVEL  ,  I  T  ,  CIFF  I 

RETURN 

ENTRY  C2FLCT 

FECTT=SNGL(H) 

WRITE  (KPRINT,5}  HEC 

WRITE     (NPRINT,6) 

WRITE     (NPRINT,7) 

RANGE( 1 J^SKGL(Co) 

RAKGE(2)  =  SNGL(CiMIN) 

RANGE(4  )  =  0 .0 

RANGE(3)=SKGL(HPLCT) 

CALL  UTPLCT  ( VE LC  ,  CEPTH  ,NPT S , RANGE  ,2 ,0 , PE CTT } 

RANGE( 1  )=i.O 

RANGE(2)=-I.O 

RETURN 

1  FCPMAT  (M',  10A8,  •   FREQUENCY  :  «,  F6  .  I  ,  •  HZ') 

2  FCRMAT  CO',  T25,  V^QDE',  14,  •    :   K  =S  Fid. 12) 

3  FCPNAT  CO',  // ,  T2j,  'NORMALIZED  t  IGENFUNCT ICN '  , 
I   •  VS  CEPTH  «) 

4  FCRMAT  COPHASE  VELGCITY  :',  F7.1,  '  M/SEC.    DSS 
^  I     l'GI57 

1  4X,  '  BOTTOM  LEAKAGE  CQEF  :',  G15.7,  //, 
$  '  GROUP  VELOCITY  :', 

2  F7.1,  '  ^'/SEC.    ITERATIONS  :',  15,  ICX, 

3  'DIFFERENCE  FUNCTION  :',  GU.E  ) 

5  FCPFAT  CI',  10A8) 

e    FORMAT  CC',  //,  T2^,  'GRAPH  OF  SCLND  VELOCITY', 

I  '  VS  DEPTH' ) 

7  FCFNAT  CO',  TI9,  'CEPTH  IN  f^ETERS.', 

1       '  SOUND  VELOCITY  IN  .METERS/SEC) 
ENC 
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SLBROUTIKE  LINEAR 

IMPLICIT  REAL*8  (A-h,  V-Z ) , 


L0G1CAL''=4  ($) 


C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<«<<<  CCMMON  eiCCKS  >>>>» 

c 


cc^^c^ 

CCMFON 
CCf^NON 

ccpyoN 

CCVMON 


CI 

D2 
D3 
C5 
D6 

CW 
II 


CCNMON  /  LL  / 


VELCdOZl),  CEPTH  (1021) 

ZZZ(  1C2I)  ,  OVZCIOCI) 

CC(5C),D£P(50 ) 

CMA)(,Cf'IN,CB,CeS(. 

H,Dh,hPLCT,UHHLF,0H2SIX,DH2C2, 

Dh2,hLGTT,CHB 

ZMAX,CSS,CECAY,CPVEL,CGVEL,GA^ 

FQV(  20)  ,FQ,  W,  V^2 

!\UMV,NMCD,NFRE(.»NNCCF,N,NP1, 

l\P21  ,KPTS,MCDE,ITTfNF2 

$GRAPF,$bGTPk,$CARLS,iCEL 


C 

C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<«<<<    FRQGRAI-     BEGINS    >»» 

C 

CNAX=C.CCQ 

CMN  =  GB 

CL=C.OCO 

CL=DEP{2) 

VL  =  CG(  I  ) 

VL=CG(2) 

CL^CU=CL-DU 

VLNVU=VL-VU 

1=2 


C 

c 


c 

c 
c 
c 


c 
c 
c 

c 


c 
c 
c 


CC    ^    J  =  I ,  i^J  P  L 
CEFJ=Dh=^-CFLCAT(J-l  ) 
DEPTH( J)=DEPJ 
IF     (DEPJ-DL)    2,2,  1 

1  l=I+i 
CL  =  DL 
VL  =  VL 

CL  =  CEP(  I  ) 
VL  =  GG(  I  ) 
CLNDU  =  L>L-DU 
VLNVU=VL-VU 

2  VELOJ=(DEPJ-DU)*VL^'VL/DLMDU+VU 

IF     (VELCJ.GT.GNAX)     CMAX=V£LCJ 
IF     (VELCJ-CMIN)     3,4,4 


■CEPTH    LCCP    BEGINS 


3    CMN  =  VELCJ 

A    VELG(J)=VELOJ 


IF     (  .NCT.SBCTFR)     GC    TO    6 
NF2=NP1+1 


•CEPTH    LCCF    ENDS 


CC    5    J=l,2C 
VELO{NPI+J)=GB 

CEPTH(NF1.+  J)  =  F+DHB*CFLGAT(  J-1) 
5  CCNTINUE 


ASSIGN  BCTTCM  VELCCITY  VALUES 


RETGRN 
ENC 
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SLfcRCUTINE  WKB 

IMPLICIT  P.EAL*8(A-h,V-Z),LQGICAL^4(S) 
C 

C       : 

C               L^ST    CH/INGED    2       7    SEPT    73         AIRY    PHASE    B.C. 
C  _ ^ 

C 

C<<<<<<<<<<<<<<<<<<«<<<<<<<<<<<<<«<<<  COFMQN  BLCCKS  >>>>>> 

C 

CCr-FCN  /  C2  /  ZZ2(iC2U,  CVZ(iOOi) 

CCf^MCN  /  D4  /  DROtCRB.DRORB 

CCNNGfN  /  L6  /  H,  CH  ,HPLCT  »  DHHLF,  Dl-2SiX,  DH2C  2  , 

1  DH2 ,hGGTT  ,DHB 

CCMMON  /  D7  /  DK,DKNIN,DK.MAX,OKCLC,DKU,DKL, 

2  DTEST,CSMAX,CE 

CCMMON    /    010/    waCB,y^CCNPl,WCCf>'IN,CVZB 

CCMiMQN    /AWKB/    DK  N  ,  D  Ii\T  ,  0  I  FF  ,  CS  ,  CN 

CCf'f^ON     /     U    /    NUFV,N^GC,NFREQ,Nf/CCF,N,NPl, 

3  NP21,KPTS,MGDE,n  »KP2 

CCMMON    /    12    /    LGKCCE,  lHNOOE,KWELL,Ii/vtLL  ,N^ELL, 
^  JHltJLGW 

CGKMCN  /  13  /  JUPPER, JBGTTt liHAP tJTPLtJTPL 

CCf^NON    /    Li    /    SGRAPF,  $6GTPR,fCAHCS,SX£L 

CCI^NON    /    PI    /    DPI,CFIC2,DPIC4,DFIRT2,D2PI 

CCMMON  /  DTH/  DTHETL 

CATA  DEX^U/l.7L2/,CWELL/6.0/ 
C 

C«<<<<<<«<<<<<<<<<<<<<<<<<<<<<<<«<<<  PROGRAM  BEGINS  »>» 
C 

CIM=O.CDC 

$TCP=. FALSE. 

JTFL=1 

CE={DKMAX-DKN)-(DKNAX+DKN) 

CKe=DS(.HT{CVZB-CE) 

IWbLL=l 

ISIG=-l 

CARG=C.CCO 

CAPGP=C.COC 
C 

C      r CEPTH  LCCP  BEGINS 

C 

c 

CC  21  J=1,NP1 

CKZ2=DE-DVZ{ J) 
C 
C     . . CFECK  FCR  CONDITIONS 

C  AND  TYPE  CF  SOLUTION 

C 

IF  (DKZ2)  1,2,3 

1  IF  (ISIG)  4,5,6 

2  IF  (ISIG)  7,21,9 

3  IF  (ISIG)  10,19,20 
C 

C      . IMAGINARY  REGION 

C 


4  CARGP=DSCRT(-DKZ2) 
CAPG=CA8G.+  CARGP 

GC  TO  21 

5  ISIG=-1 
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GC  TO  ^ 


C 

c 
c 


c 
c 
c 


IC 


11 


13 
1^ 


15 


16 


17 
IS 


19 


2C 


C 
C 
C 
C 


CyiRGP  =  CSGRT(-CK22J 

CX=DKZ/ (CKZ+DAPGP) 

C1M=LIM  +  DKZ*(DX-1  .OCO  J'i^O  .  5D0*Ch 

ISIG=-L 

JTFL=J-1 

C/RG=(2.0CO-CX)^D/iRGP*C.5DO 

GC    TO    21 


■ZERO    VERTICAL    WAVE    NUMBER 


7     I£ACK=-1 
6    ISIG=G 
GC    TO    2  1 

S    IEACK=1 
GC    TO    8 


■REAL    REGICN 


IS 
IF 
JT 
CK 
CX 
C/ 
IF 
C2 
CT 
CI 
$T 
GC 
CT 
IF 
IF 
CA 
IVs 
JT 
CT 
CI 
GC 
A  = 
E  = 
CE 
CE 
AC 
EC 
CI 
IF 
CI 
GC 
CI 
CI 
CA 
GC 


IG  =  1 

{5T0P)    GO    TO     11 
FU=J 


Z  =  DSQ 
=  CARG 
RG=CA 

($TQ 
=  2.CC 
hETU  = 
NT=DT 
CP=.T 

TO  1 
1  =  DMC 

(DAR 

(I  WE 
PG=C. 
ELL=I 
Flj=J 
FETU= 
NT  =  DT 

TO  1 
CSIN{ 
CCCS( 
L  =  CEX 
fL=CE 
=  DEL  = 
=  DEL-^ 
NTP=D 

(DIN 
NTP=D 

TO  1 
NT  =  DI 
NT=DI 
RG=0. 

TO  2 


ST(CKZ2) 

P/(CKZ+DARGP) 

F.G^CAFGP'i^(DX-1.0C0}*0.5DC 

Pi     GO    TG    13 

C-D£XP(2.GDC-CARG^DH  ) 

CAT  AN ( (D2-1  .CD0)/(D2+1.0CC)) 

F-tTC 

RUE. 


C(DINT,D2P1  ) 
G-'Ch-CWELL  ) 
LL.EG.KWELL) 

ccc 

WELL+1 


15,  1^,  lA 

GC    TO    28 


CPIO- 
FETU 

e 

CTl) 
CTl) 

F(DAI 

XP(-I 

(A+E 

(A  +  e 

INT-DTl  +  DAT/iN2(AC,B0) 

TF.GE.CINT-CTEST)    GC    TO 

INTP- 

6 

NTP 

NT+12.0D0-CX)*DKZ*C.5D0*DH 

ceo 
1 


iRG=^DH) 

■DARG-^Ch  ) 
)  +  CEiML"(A-B) 
)-DEML-(A-B) 

•DTl  +  DAT/iN2(ACt 
E.CINT-CTEST) 

'+0PI02 


17 


CKZ  =  DS(.RT(CKZ2) 
ISIG=+1 

CINT=DINT.+  DKZ=*OH 
IF     ($TCF)     GO    TO    13 
JTPU=J 
GC    TO     12 

CKZ=CSQRT(DKZ2) 
CINT=DINT+CKZ*DH 


21    CCNTINUE 


•DEPTH    LCCF    ENDS 
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c 

c 
c 


c 
c 


IF  (DKZ2)  22,27,25 


■DECAYING  SCLtTICN 


12    CKS=DARGP 

C/iFG=(C/5p,G-C.500*DARGP  )^DH«2.0DC 

IF  (CARG-DEXPU)  2b,24,2A 
2  3  C5=(DKS-DRCRb*^DKBl*0EXF(-DARG)-«2.0CO/(OKS+CPCPE=PDKe) 

DThETL=CATAN(( l.0D0+D3)/{  1.0D0-D2)) 

GC  TO  27 
2^  CTFETL=DPIQ4 

GC  TO  27 


JTFL=NFI 

IF  (DKB.FC.O.ODG)  GC  TO  26 
CTFETL=Ca*TAN(DKZ/  (CK8*CRQRB)  ) 
GC  TO  27 


ECTTOM  REFLECTED 


It    CTf-ETL=CFIC2 


COMPUTE  DIFFERENCE 


26 


CS=CINT+DTFETL 
CINT=DINT-DThETU 
DIFF  =  DS-CF'>DPI 
MnELL  =  I  i\ELL 
RETURN 

Il>ELL  =  IyNElL+l 
GC    TO    2^ 

END 
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